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tion  process  is  implemented  by  using  multigrid  V-cycles;  the  iterative  and  multilevel 
results  are  cotnpared  and  discussed  in  detail.  The  computational  savings  resulting 
from  the  multilevel  process  are  then  discussed. 


DISCLAIMER 


The  computer  program  in  Appendix  B  is  supplied  on  an  “as  is”  basis,  with  no 
warrantees  of  any  kind.  The  author  bears  no  responsibility  for  any  consequences  of 
using  this  program. 


Vll 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

II.  PROBLEM  STATEMENT .  3 

A.  BACKGROUND .  3 

B.  POINT-SOURCE  PROBLEM .  5 

III.  DISCRETIZATION .  7 

A.  BASIS  FUNCTIONS:  FINITE  ELEMENTS .  12 

B.  EVE  STENCILS  .  15 

1.  Example;  2— D  Potential  Flow .  15 

2.  .Application  to  Axisymmetric  Heat  Equation .  20 

IV.  RELAXATION  SCHEMES .  27 

A.  ITERATION  MATRICES .  28 

B.  ALTERNATIVE  SCHEMES:  LINE  RELAXATION .  35 

C.  LOCAL  MODE  ANALYSIS .  38 

V.  ITERATIVE  SOLUTION . .  53 

VI.  MULTILEVEL  SOLUTION .  61 

A.  INTERGRID  TRANSFERS .  67 

B.  AMPLIFICATION  MATRIX  .  73 

C.  NUMERICAL  RESULTS .  83 

VII.  CONCLUSION .  95 

APPENDIX  A.  THE  DISCRETE  FOURIER  TRANSFORM  ....  99 

APPENDIX  B.  INTERPOLATION  TOOLS . 105 

REFERENCES . 109 

INITIAL  DISTRIBUTION  LIST  . Ill 


IX 


I.  INTRODUCTION 


A  topic  of  general  interest  is  the  numerical  solution  of  partial  differential  equa¬ 
tions  (PDEs),  such  as  the  Navier-Stokes  equation  (Equation  II. 1).  Many  aspects  of 
tl)e.se  types  of  problems  present  interesting  challenges  to  constructing  numerical  so¬ 
lutions,  such  as  non-linearities,  high-speed  flows  that  cause  numerical  peculiarities, 
unusual  boundary  conditions,  uncommon  geometries,  and  so  on.  Developing  pro¬ 
cesses  that  solve  these  types  of  problems  is  difficult,  and  often  results  in  complicated 
solutions  schemes  which  are  costly  to  implement.  A  common  goal  is  to  attempt  to 
streamline  existing  solution  processes,  or  develop  new  ones,  in  order  to  reduce  com- 
])utational  complexity  and  cost. 

At  issue  is  how  to  transform  a  continuum  problem  into  a  discrete  one  (dis¬ 
cretization),  and  how  to  construct  a  numerical  process  that  will  solve  the  discrete 
problem.  The  finite  volume  element  (EVE)  method  (see  [Ref.  1])  has  proven  to  be 
a  useful  tool  in  developing  discretizations,  particularly  for  problems  that  require  the 
enforcement  of  conservation  laws.  One  of  the  more  successful  methods  for  streamlin¬ 
ing  the  solution  of  PDEs,  particularly  elliptic  equations,  is  a  multilevel  technique,  to 
wit,  inultigrid  (see  [Ref.  2],  [Ref.  3],  [Ref.  1]).  While  this  method  enjoys  remarkable 
success  in  solving  certain  classes  of  both  linear  and  non-linear  PDEs,  its  applicability 
to  solving  other  types  of  equations  awaits  development.  The  goal  of  this  work  is  two¬ 
fold:  to  apply  the  EVE  method  to  discretize  a  particular  equation,  and  to  implement 
multigrid  in  solving  the  resulting  discretized  problem. 

The  approach  that  is  taken  is  to  begin  with  a  specific  physical  problem  which 
results  in  the  Navier-Stokes  equation,  and  consider  a  subsidiary  problem,  namely  the 
heat  equation;  due  to  the  nature  of  the  problem,  cylindrical  coordinates  are  used. 
The  resulting  problem,  to  which  the  EVE  discretization  and  multilevel  solution  are 
applied,  is  the  axisymmetric  heat  equation,  involving  the  use  of  a  point  source.  More 
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s])ecifically,  finite  differences  are  used  to  discretize  the  temporal  portion  of  the  prob¬ 
lem.  resulting  in  the  use  of  a  fully  implicit  backward  time-stepping  scheme,  and  the 
FVE  method  is  used  to  discretize  the  spatial  part.  The  application  of  the  FVE 
method  to  a  problem  in  cylindrical  coordinates  is  new,  and  results  in  stencils  which 
are  analyzed  extensively.  Once  the  discretization  has  been  constructed,  a  number  of 
relaxation  schemes,  including  both  .Jacobi  and  Gauss-Seidel,  are  considered;  a  thor¬ 
ough  analysis  of  these  schemes  is  done,  using  both  the  spectral  radii  of  the  iteration 
matrices  and  local  mode  analysis.  The  goal  is  to  choose  one  of  the  relaxation  schemes 
for  use  in  an  iterative  solution  of  the  problem;  Gauss-Seidel  is  the  method  of  choice. 
The  specifics  of  the  multilevel  technique  are  then  developed  and  analyzed.  Local 
mode  analysis  is  performed  on  the  components  of  the  amplification  matrix,  resulting 
in  the  two-level  convergence  factors  for  various  combinations  of  the  operators.  The 
multilevel  solution  process  is  implemented  by  using  multigrid  V-cycles;  the  iterative 
and  multilevel  results  are  compared  and  discussed  in  detail.  The  computational  sav¬ 
ings  resulting  from  the  multilevel  process  are  then  discussed.  Insofar  as  multigrid  is 
cpaite  successful  in  a  number  of  different  contexts,  the  results  of  our  work  are  some¬ 
what  disappointing  in  that  we  are  not  yet  able  to  achieve  the  same  level  of  success. 
On  a  more  positive  note,  this  work  applies  the  FVE  method  to  a  problem  in  cylindri¬ 
cal  coordinates,  and  makes  use  of  Maple  software  to  compute  the  necessary  integrals 
to  construct  the  discretization  (see  Chapter  III  and  Appendix  B).  Additionally,  the 
techniques  that  are  applied  do  result  in  a  solution  to  the  problem. 
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II. 


PROBLEM  STATEMENT 


A. 


BACKGROUND 

Our  ultimate  goal  is  to  be  able  to  numerically  solve  the  Navier-Stokes  equation. 


-T — h  u  ■  Vu  =  — Vp  + 
at  p 


(II.1) 


where  u  is  the  velocity  vector,  t  is  the  time,  p  is  the  density,  p  is  the  pressure,  and  u 
is  the  kinematic  viscosity.  One  approach  to  such  solutions  is  to  consider  this  problem 
as  a  collection  of  smaller  problems.  By  solving  each  of  the  smaller  problems,  we 
ran  assemble  the  collection  of  solution  processes  to  solve  the  original  problem.  The 
|)urpose  of  this  work  is  to  focus  on  an  initial  step  in  the  eventual  construction  of  a 
numerical  solution  to  Equation  II.  1. 

One  way  to  subdivide  a  problem  such  as  this  into  smaller  problems  is  to 
consider  a  specific  physical  problerri  that  gives  rise  to  the  Navier-Stokes  equation,  and 
then  isolate  the  different  aspects  of  that  problem.  To  that  end,  we  consider  a  semi¬ 
infinite  block  of  metal,  with  a  heat  source  applied  to  the  horizontal  surface;  above  the 
liorizontal  surface  is  an  (inviscid)  nonconducting  gas.  The  result  is  a  pool  of  molten 
metal  with  a  free  surface,  surrounded  by  the  solid  portion  of  the  unmelted  block  (see 
Figure  1).  The  total  heat  flux  Q  is  constant,  and  far  away  the  solid  approaches  the 
uniform  cold  temperature,  Tc.  The  resulting  thermal  and  flow  fields  are  assumed  to 
be  axisymmetric  and  steady,  and  are  governed  by  conservation  of  energy  in  the  solid 
and  l>y  conservation  of  energy,  momentum,  and  mass  in  the  pool.  We  therefore  have 
the  following  system  of  equations: 


solid 

liquid 


ri'T' 

—  -b  u  •  VT  = 

at 

du  ^  1 

at  p 

V  •  u  =  0 


(11.2) 

(11.3) 

(11.4) 

(11.5) 
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where  T  is  the  temperature,  k  is  the  thermal  diffusivity,  and  u,  t,  p,  p,  and  u  are  as 
above. 


Q 


solid 


Figure  1.  The  molten-pool  problem. 

The  portion  of  the  above  problem  which  is  our  focus  is  the  conduction  of  heat 
in  the  solid,  governed  by  the  heat  equation  (Equation  II.2),  to  which  the  following 
boundary  conditions  apply: 


surface  z  =  0  : 

dT 

(11.6) 

axis  r  =  0  : 

o 

II 

5  1  ^ 

(11.7) 

far  away  z  — >•  oo  : 

: 

(11.8) 

where  k  is  the  thermal  conductivity,  and  q{r)  is  the  imposed  surface  heat  flux  (large 
at  r  =  0,  falling  off  to  zero  at  some  small  value  of  r,  such  that  q[r)2'Krdr  =  Q);  the 
condition  in  Equation  II. 7  results  from  the  axisymmetric  nature  of  the  problem.  Since 
the  ultimate  goal  is  to  solve  the  Navier-Stokes  equation  in  cylindrical  coordinates,  the 
heat  equation  will  also  be  solved  in  this  coordinate  system. (see  [Ref.  4]) 
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B.  POINT-SOURCE  PROBLEM 


Our  goal  is  to  take  the  first  step  toward  solving  the  Navier-Stokes  equation 
(Equation  II. 1)  that  arises  from  considering  the  molten-pool  problem  by  solving  the 
associated  heat  equation  (Equation  II. 2).  We  therefore  assume  cylindrical  geome¬ 
try  and  azimuthal  symmetry.  We  further  simplify  the  conditions  imposed  on  Equa¬ 
tion  1 1.2  by  considering  that  the  heat  source  applied  to  the  horizontal  surface  of  the 
block  is  a  point  source,  with  the  total  heat  flux  (5  =  1,  and  that  far  away  the  solid 
ai)i)roaches  the  uniform  cold  temperature  of  Tc  =  0.  That  is,  the  imposed  heat  flux 
(I[t)  =  S{r)  is  the  Dirac  delta  function,  with  q(r)27vrdr  =  1. 

The  existence  of  the  solution  to  this  type  of  problem  is  well  established.  How¬ 
ever.  as  is  the  case  in  general  when  Neumann  boundary  conditions  apply,  absent  a 
specified  initial  temperature  distribution,  the  solution  is  not  unique  (see  [Ref.  5]). 
While  we  are  interested  in  solving  the  point-source  problem  numerically  in  cylindrical 
coordinates,  it  has  spherical  symmetry,  so  it  can  be  solved  analytically  in  spherical 
coordinates  more  easily.  Therefore,  for  the  analytic  solution,  we  rewrite  Equation  II.2 
in  spherically  symmetric  coordinates: 


dt  R2dR\  dRj' 


(11.9) 
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where  erfc(x)  is  the  complementary  error  function, 


erfc(x)  =  1  —  erf(x),  erf(x)  =  —7=  [  e 

x/TT  Jo 


2  2 

-7=  /  e-"  ds. 


As  t  — )■  00,  the  steady  solution  becomes 


T{R) 


27TkR 


(11.14) 


(for  a  derivation  of  this  result,  see  [Ref.  5]).  In  (r,  z)  coordinates,  where  R  =  vrM-?, 
these  become 


/j.2  I  2)2 


2v^  }  ’ 


r(r,z)  = 


(ir.i5) 

(11.16) 


For  computational  purposes,  the  idealized  problem  of  a  semi-infinite  solid 
is  truncated  to  a  finite  domain  in  cylindrical  coordinates  with  azimuthal  symme¬ 
try.  .At  the  far  boundaries  on  this  finite  domain,  homogeneous  Dirichlet  conditions 
are  imposed  to  approximate  the  conditions  in  the  unbounded  problem.  To  further 
simplify  computational  work,  we  assume  that  the  {r,z)  domain  is  the  unit  grid 
=  [0,  l]x[0, 1]  €  7^'^;we  will  eventually  assume  k  =  1.  Thus,  the  equation  we 


wish  to  solve  is 


kV^T, 


(11.17) 


with  boundary  conditions 


OT 

surface  z  =  0  :  =  —S(r) 

az 

dT  , 

axis  r  =  0  :  =  0 

or 

far  boundaries  r.z  =  I  :  T  =  0. 


(11.18) 

(11.19) 

(11.20) 
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III. 


DISCRETIZATION 


In  order  to  solve  the  conduction  problem  numerically  (Equation  II. 2,  with 
l)ouudary  conditions  Equations  11.18,  11.19,  and  11.20),  the  equation  representing  the 
continuum  problem  must  be  discretized.  That  is,  the  values  of  the  unknown  function 
T  are  to  be  determined  from  a  set  of  equations  which  in  some  sense  approximate 
Equation  II. 2.  Ultimately,  a  discrete  approximation  to  the  heat  equation  should  meet 
the  following  criteria: 

1)  F’rovide  for  a  unique  solution  to  the  problem. 

2)  The  solution  should  be  “close”  to  the  exact  solution  for  “sufficiently  small” 
grid  spacings. 

.’1)  The  solution  should  be  “effectively  computable.” 

The  significance  of  property  1)  is  obvious;  property  2)  relates  to  the  questions  of 
convergence  and  consistency  for  the  discretization  scheme.  Property  3)  relates  to: 
a)  the  amount  of  computational  work  required  to  solve  the  problem,  and  b)  the 
behavior  of  roundoff  errors  in  the  computed  solution.  The  growth  of  the  roundoff 
ei  ror  is  related  to  the  notion  of  the  stability  of  the  discretization  scheme,  (see  [Ref. 

fi]) 

Solution  of  the  heat  equation  (Equation  11.2)  requires  treatment  of  both  space 
and  time  derivatives.  The  Finite  Volume  Element  (EVE)  method  is  used  to  discretize 
the  spatial  portion  of  the  problem;  finite  differences  are  used  to  discretize  the  tem¬ 
poral  portion.  The  approximation  properties  of  the  EVE  method  are  fundamentally 
different  from  those  associated  with  the  finite  difference  method.  Hence,  this  ap- 
|uoach  treats  time  and  space  differently.  The  goal  in  applying  the  EVE  method  to 
produce  spatial  discretization  is  to  make  use  of  the  advantages  of  the  Finite  Volume 
(FV)  method,  its  ability  to  be  faithful  to  the  physics  in  general  and  conservation 
in  particular,  coupled  with  the  flexibility  of  the  flnite  element  representation  of  the 
unknown  functions  (see  [Ref.  1]).  As  outlined  in  [Ref.  1],  the  basic  approach  is 
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1,0  [partition  the  spatial  grid  in  two  ways:  as  the  union  of  a  set  of  finite  elements, 
the  vertices  of  which  comprise  the  grid  points  on  which  the  unknowns  are  defined 
(see  Figure  2),  and  as  the  union  of  a  set  of  control  volumes  (see  Figure  3),  one  for 
each  grid  point  (the  boundary  points  require  separate  treatment,  as  will  be  discussed 
later).  The  elements  are  used  to  form  basis  functions,  a  linear  combination  of  which  is 
used  to  approximate  the  unknown  function.  Upon  substitution  of  the  approximation 
into  the  equation  to  be  solved,  integrals  over  each  control  volume  are  taken.  The 
integrals  enforce  conservation  on  each  control  volume,  and  therefore  the  partition 
enforces  conservation  on  the  entire  domain.  The  system  of  equations  generated  by 
integration  yields  the  discretization  of  the  problem. 


0  1  2  •••  N 

r  _ ^ 

Figure  2.  Partitioning  of  problem  domain  cross-section  fi  by  elements. 

The  process  of  discretizing  Equation  II. 2  in  both  space  and  time  is  complicated, 
the  application  of  the  FVE  method  to  the  space  derivatives  being  the  more  complex 
portion  of  this  process.  Therefore,  before  proceeding  to  a  discussion  of  the  FVE 
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Figure  3.  Partitioning  of  problem  domain  cross-section  0,  into  sub-volumes. 


method,  the  finite  difference  method  is  applied  to  the  time  derivative.  In  order  to 
facilitate  this  presentation,  we  present  a  one-dimensional  example  of  this  type  of 
discretization  (see  [Ref.  7]).  Consider  the  one-dimensional  version  of  Equation  II. 2, 


dT  _  d^T 

'm  ~ 


(III.l) 


where  k  =  1.  One  finite-difference  approximation  to  Equation  III.l  is 

2fc,n+l  ~  Tk,n  _  7fc+l,n  —  2Tfc_„  -f  Tk-\,n 
~g  “  P  ’ 

where  T  is  the  exact  solution  of  the  approximating  difference  equations. 


Xk  —  kh  &nd  tn  —  gn,  (fc,  n  =  0, 1, 2,  •  •  •). 


Let  r  —  and  this  can  be  written  as 


Tk,n+\  —  '>'Tk-l,n  ~  (1  —  '^1')Tk,n  +  'l'Tk-\-\,n 
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Thus,  we  can  compute  the  unknown  Tk^n+i  at  the  (n  +  l)st  time  step  using  the 
known  values  from  the  nth  time  step,  giving  an  explicit  formula  for  determing  the 
unknowns. [Ref.  7] 

While  this  explicit  relation  provides  a  simple  means  to  compute  unknown 
values,  is  has  a  serious  drawback.  The  process  is  only  valid  for  0  <  <  5.  or 

•j<'i  ,  and  therefore  the  time  step  =  g  \s  necessarily  small.  Moreover,  in  order 
to  achieve  reasonable  accuracy,  Ax  =  h  must  also  be  kept  small.  There  is,  however, 
a  method  which  reduces  the  total  amount  of  calculation  and  is  valid  (i.e.,  convergent 
ajid  stable)  for  all  finite  values  of  i\  which  was  proposed  by  Crank  and  Nicholson 
(see  [Ref.  7]).  The  technique  is  to  consider  Equation  III.l  as  being  satisfied  at 
the  midpoint  {A,7r,  (n  +  ^)g}  and  replace  by  the  average  of  its  finite-difference 
approximations  at  the  nth  and  (n  -f  l)st  time  steps.  That  is,  the  equation 


is  approximated  by 


Tk-l,n+l  —  ‘^Tk,n+l  +  Tk-l,n  —  ‘2Tk,n  + 


which  gives 


—  'l'Tk-\,n+\  -k  (2  -f  2r)Tk,n+l  ~  f'Tk+X^n+l  —  f’Tk-l,n  +  (2  —  2r)Tk,n  +  1'Tk+l,ni 

w  her  e  —  11^2  •  Therefore,  we  now  have  three  unknowns  at  the  (n  -f-  l)st  time  step 
written  in  terms  of  three  known  values  at  the  rrth  time  step.  Thus,  the  Crank- 
Nicholson  method  generates  a  set  of  equations  which  must  be  solved  simultaneously; 
it  is  an  implicit  method. 

Although  the  Crank-Nicholson  method  for  Equation  III.l  is  stable  for  all  pos¬ 
itive  values  of  r  in  the  sense  that  all  errors  eventually  tend  to  zero  as  n  tends  to 
infinity,  large  values  of  r,  e.g.,  40,  can  introduce  oscillations  in  the  solution  (see  [Ref. 
7]).  Therefore,  a  more  general  finite-difference  approximation  to  Equation  III.l  is 


found  by  using  a  weighted  average  of  the  finite-difference  approximations  of  at 
tlie  ?i-tli  and  (n  +  l)st  time  steps,  which  is  given  by 


-  2r,. 


-t-(l-a) 


'2Tk,n  +  Tk+i^n 


where  0  <  a  <  1  is  possible.  Note  that  a  =  0  gives  the  explicit  scheme,  a  =  Crank- 
Nicholson,  and  a  =  1  a  fully  implicit  backward  time-difference  method.  The  scheme 
is  unconditionally  valid  (stable  and  convergent)  for  |  <  a  <  1,  but  for  0  <  «  <  |  we 
must  have  7-  <  1/(2  —  da).  Thus,  for  example,  for  cv  =  0,  r  =  ^  <  d,  or  f/  <  '>s  a 

requirement  to  guarantee  validity. [Ref.  7] 

We  now  apply  this  finite-difference  method  of  discretizing  the  time  derivative 
to  Equation  11.2.  Following  the  work  done  above,  a  general  finite-difference  approxi¬ 
mation  for  the  time  derivative  of 


dT 

It 


kV^T 


(in.2) 


is  given  by 


=  qkV*T”+‘  +  (1  -  a)KV’r”, 


:iii.3) 


where  g  is  the  time  step,  and  0  <  cr  <  1  allows  for  a  weighted  average  of  the  current 
and  future  time  steps.  The  current  time  step  is  designated  by  n,  ?^  +  1  is  the  next 


time  step;  the  value  of  a  determines  the  type  of  time  stepping.  However,  we  have 
determined  that  in  order  to  guarantee  the  validity  (stability  and  convergence)  of  the 
explicit  version  of  this  scheme  (a  =  0),  we  must  have  g  <  y-  consider  this 
value  of  g  to  be  too  small  to  be  of  computational  use  (which  is  verified  by  numerical 
experiments),  and  therefore  we  will  use  the  values  o:  =  ^  and  a  =  \.  Thus,  we  have 
the  semi-discrete  relationship  between  the  current  and  subsequent  time  steps: 


(1  _  crKV2)r"+i  =  (i  -I-  (1  _  a)«V2)T".  (III.4) 
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This  relationship  is  used  to  establish  a  discrete  set  of  equations,  which  gives  rise  to 
the  discretization  of  the  problem. 


A.  BASIS  FUNCTIONS:  FINITE  ELEMENTS 

We  begin  the  discussion  of  discretizing  Equation  II. 2  in  the  spatial  variables  by 
recognizing  that  it  is  necessary  to  make  an  assumption  about  what  types  of  functions 
may  be  used  to  approximate  the  unknown  function  T;  we  consider  that  T  can  be 
ajjproximated  by  continuous,  piecewise  linear  functions.  A  basis  is  then  chosen  for 
the  space  in  which  these  functions  are  found.  That  is,  a  set  of  functions  is  selected 
whose  linear  combinations  determine  the  space  of  functions  of  interest.  Therefore,  we 
choose  basis  functions,  (pkji’’’,  which  are  assumed  to  be  continuous  and  piecewise 
linear,  with 

^■=o  i=o 

The  next  step  is  to  determine  how  to  construct  these  basis  functions;  the  element 
partition  of  the  domain  Q,  is  used  as  a  framework  for  this  purpose,  (see  Figure  2). 

Since  we  are  assuming  cylindrical  geometry  with  azimuthal  symmetry,  the 
directions  of  interest  are  r  and  z\  a  uniform  grid  of  step  size  h  =  {N  =  number 
of  grid  points)  is  applied  in  both  directions.  In  this  way,  we  can  designate  function 
values  on  the  grid  17  by 

Tk,j  =  T{rj,Zk), 

where  Zk  —  kh,  Vj  =  j/Ti,  and  k,j  are  the  row  and  column  indices.  Each  square 
is  subdivided  into  two  triangular  elements  by  a  diagonal  (oriented  in  the  direction 
of  increasing  v  +  zY;  the  basis  functions  to  be  used  in  approximating  the  unknown 
function  are  constructed  using  these  triangular  elements.  For  each  interior  grid  point, 
or  node^,  there  are  six  associated  triangular  elements:  NNE,  ENE,  SE,  SSW,WSW, 
NW  (see  Figure  4). 

'The  orientation  of  the  elements  is  based  on  the  necessity  to  apply  this  technique  to  solving  the 
rnolten-pool  problem.  In  particular,  we  will  be  required  to  track  the  movement  of  the  phase-change 
boundary  between  solid  and  liquid;  element  orientation  has  been  chosen  to  facilitate  keeping  track 
of  this  moving  boundary. 

^For  a  two-dimensional  grid  with  N  intervals  in  each  dimension,  there  are  (N  —  1)^  interior  grid 
points,  {N  -t-  1)^  total  nodes. 
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(k+lj+I) 


(k+l,j) 


Ck-lJ-l)  (k-lo) 

e  (p  ^  ^ 

Figure  4.  Hat  function,  top  view,  k  is  the  row  index  [z  direction),  j  is  the  column 
index  {r  direction),  direction  is  into  the  plane  of  the  paper  (©). 

If  we  suppose  that  this  collection  of  elements  is  joined  along  shared  borders, 
and  that  the  elements  can  be  ‘^stretched”,  then  the  finite  element  for  a  particular 
node,  say  (A;,  j),  is  formed  by  anchoring  the  collection  along  its  external  boundary  (in 
the  (r,  z)  plane)  and  extruding  the  grid  point  in  a  direction  perpendicular  to  the  (r,  z) 
plane  and  parallel  to  the  tangent  to  the  ip  direction  at  the  grid  point  (A;,  j),  forming 
the  familiar  ‘4aat”  function  (see  Figure  5).  Note  that 

f 

1  at  the  node  (r^, 

0  at  all  other  nodes 

and  is  piecewise  linear  over  each  triangular  element,  giving  the  hat  shape. 

By  choosing  these  hat  functions  to  be  our  basis  functions,  the  coefficients  of 
the  z)  in  Equation  IIL5  become  the  values  I3kj  =  T/tj.  That  is,  we  can  now 

specify  the  nodal  values  of  T  on  the  grid  as  Tkj  =  T{rj,Zk),  with 

T(r,z)  =  ^^Ti,jMr,z).  (111.6) 

k=zO  j=0 
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Figure  5.  Hat  function,  oblique  view. 

Since  the  are  piecewise  linear,  a  continuous  representation  for  function  values  of 
T  is  found  by  linear  interpolation  of  nodal  values  between  nodes.  For  example. 


Tk,j  + 


Up  to  now,  we  have  considered  the  sampled  values  of  the  unknown  T  at  the 
{N  +  1)'^  grid  points  as  an  {N  +  l)x(A^  +  1)  two-dimensional  array;  it  will  later 
l)e  convenient  to  consider  T  as  a  one  dimensional  array.  One  simple  method  to 
accomplish  this  is  to  order  the  elements  Tkj  lexicographically.  For  example. 


7o,o  ^0,1 

7\,o  Ti^i 


To,n+i 

Ti,n+i 


can  be  written  as 


Tjyj+ifi  7V+i,i  •  •  •  2V+i,/v+i 


(iro,o,  ?0,i,  .  .  .  ,  To,Ar+l,  Ti^o,  Ti^i,  .  .  .  ,  ri,;v+l5  •  •  •  ?  TW+l.O)  TV+i,!, . . . ,  TN+i,N+iy 
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vvliich  we  will  normally  consider  to  be  a  column  vector.  Thus,  we  can  relabel  the 
values  Tfcj  =  Ti,  I  —  {N  +  l)k  +  j  +  I  and  A;,  j  =  0  :  +  1,  in  Equation  III. 6  so  that 

{N+if 

r(.',-)=  r, *(>•,»).  (in.7) 

/  =  1 

B.  FVE  STENCILS 

With  the  construction  of  the  basis  functions  over  the  finite  elements,  we  can 
now  incorporate  the  control  volumes  to  complete  the  finite  volume  element  discretiza¬ 
tion.  Below  is  an  example  the  FVE  method  applied  to  a  simpler  problem;  the  purpose 
is  to  motivate  and  explain  the  work  that  is  necessary  to  handle  the  specifics  of  apply¬ 
ing  the  FVE  method  to  Equation  II. 2.  This  explanation  closely  follows  the  example 
in  [Ref.  1]. 

1.  Example:  2— D  Potential  Flow 

The  example  problem  we  consider  is  the  potential  flow  problem  in  Euclidean 
two-dimensional  coordinates.  The  domain  for  the  problem  is  the  unit  grid  = 
[0,  l]x[0, 1]  G  71^1  where  x  =  0  and  y  =  0  are  the  near  boundaries  and  a;  =  1  and 
y  =  1  are  the  far  boundaries.  The  governing  equation  is 

V .  yvi,)  = ,,  (in.8) 

where  p  is  the  density,  ip  is  the  flow  potential  (i.e.,  ipx  and  ipy  are  the  components 
of  the  fluid  velocity  in  the  x  and  y  directions),  and  rj  is  the  interior  source  flow  rate. 
The  boundary  conditions  are 

near  boundaries  :  [pVip)  ■  h  =  ipi  (Neumann  boundary  condition), 
far  boundaries  :  ip  =  ipQ  (Dirichlet  boundary  condition), 

where  n  is  the  outward  normal.  Suppose  that  V  is  one  of  the  control  volumes  in  fl, 
similar  to  those  depicted  in  Figure  3,  where  r  and  2  are  replaced  by  x  and  y.  Since 
this  problem  is  in  two  dimensions,  the  control  volumes  V  are  actually  areas,  and  the 
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corresponding  surfaces  S  are  the  perimeters  of  these  areas.  Integrating  Equation  III. 8 
ovei'  V  ,  we  have 

f  V-ipW'iP)dV=  /  pdv. 

Jv  JV 

By  applying  the  Gauss  Divergence  Theorem,  the  volume  integral  on  the  left-hand 
side  becomes  a  surface  integral, 

^(pVv)  •  hdS  =  ^  ridV.  (III.9) 

In  terms  of  the  physical  units  associated  with  the  potential  flow,  the  relationship  is 

mass  length  mass 

- — ^ - -  area  =  - - - -  volume; 

volume  time  time  •  volume 

each  side  represents  a  flow  rate  of  mass  per  unit  time,  and  {pVi/))  ■  h  represents  a  flux 
across  S.  Therefore,  the  net  flow  from  the  interior  source  in  V  is  balanced  by  the 
flow  across  the  surface  S,  which  can  be  interpreted  as  a  mass  conservation  law  for  the 
control  volume  E.fRef.  1] 

By  partitioning  the  domain  D  into  a  finite  collection  of  control  volumes  (see 
Figure  3,  where  r  and  z  are  replaced  by  x  and  y),  we  can  impose  on  each  control 
volume  the  condition  represented  in  Equation  III. 9  (the  boundaries  require  special 
treatment;  a  discussion  of  these  considerations  is  included  later).  As  noted  earlier, 
the  imposition  of  conservation  on  each  control  volume  results  in  conservation  over 
the  entire  domain.  This  integral  condition  on  each  control  volume  will  be  used  to 
compute  the  discretization;  the  number  of  discrete  equations  is  the  number  of  control 
volumes,  say  m,  used  to  partition  D.  To  complete  the  discretization  process,  the 
unknown  function  0  will  be  approximated  by  (d  in  7?.™.  By  appropriately  choosing 
basis  functions,  we  can  construct  an  approximation  v  expressed  in  terms  of  its  nodal 
values.  Consider  a  triangular  element  partition  similar  to  that  in  Figure  2  (where  r 
and  0  are  replaced  by  x  and  y)  and  let  T  be  the  space  of  continuous,  piecewise  linear 
functions  on  f2  associated  with  this  partition.  Ignoring  for  the  moment  conditions  at 
the  boundaries,  the  EVE  discretization  of  Equation  III. 8  comes  from  finding  a  u  G  T 
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siicli  that  Equation  IIL9  holds  for  each  of  the  control  volumes  in  the  partition  of 
The  problem  in  is  then  defined  when  v  is  expressed  in  terms  of  nodal  basis  for  T: 

m 

u  =  ^  (III.IO) 

u;=l 

where  is  the  value  of  v  at  the  roth  node  and  (j)^,  is  the  “hat  function”  associated  with 
the  u’th  node  (as  above,  we  lexicographically  order  the  nodes,  resulting  in  a  single 
subscript).  Substituting  Equation  III.IO  into  Equation  III. 9,  we  have  the  matrix 
('(|ualion 

Lu  =  f,  (III.ll) 

where  L  is  the  nxn  matrix  with  entries 

Uz=  i  {pV4>.)-ndS  (III.12) 

and 

[  rjdV  (III.13) 

(except  for  boundary  conditions). [Ref.  1] 

The  treatment  of  the  boundaries  depends  the  type  of  boundary  condition. 
Neumann  conditions  can  be  imposed  indirectly  by  substituting  the  flux  value  V’l  into 
the  appropriate  term  in  Equation  III. 9.  Specifically,  suppose  we  choose  the  control 
volume  V  that  includes  the  origin  (lower  left-hand  corner.  Figure  3),  where  the  surface 
of  the  control  volume  coincides  with  the  boundary  along  the  two  axes.  The  integral 
condition  for  V  is 

f  il^idS  -h  /  (pVV’)  -ndS^  [  TjdV-,  (III.14) 

Jx,y—0  Jx^y:^0  Jv 

Equation  III.  14  is  substituted  for  the  interior  condition  in  Equation  III. 9  into  the 
discrete  approximation  v  for  the  corner  volume  containing  the  origin.  [Ref.  1] 

Dirichlet  conditions  are  imposed  directly;  Uyj  takes  on  the  value  of  ■01  at  each 
Dirichlet  node,  i.e.,  each  node  on  the  far  boundaries  {x,y  =  1).  This  approach  may 
lead  to  fewer  unknowns  than  equations,  a  problem  easily  resolved  by  discarding  the 
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(-'<|uatioiis  associated  with  the  Dirichlet  nodes.  In  this  case,  the  collection  of  control 
volumes  no  longer  partitions  the  domain  and,  therefore,  conservation  no  longer  is 
strictly  enforced.  However,  it  is  at  the  Dirichlet  nodes  where  the  loss  of  conservation 
is  not  a  concern  in  general. [Ref.  1] 

Assuming  that  there  is  a  balance  in  the  number  of  equations  and  unknowns,  we 
now  must  implement  a  numerical  rule  for  evaluating  the  integrals  in  Equations  III.  12 
and  111.13.  For  Equation  111.12,  we  use  the  following  rule  on  each  segment  A  that  is 
part  of  the  interior  surface  S  associated  with  a  volume  V: 


ipVS)  ■  hdS  cii  p[P a)  [  V</>  •  ndA 

1  JA 


where  Pa  is  point  of  intersection  of  A  and  the  grid  lines  passing  through  the  node 
associated  with  V .  For  Neumann  boundary  segments,  we  use  the  quadrature  rule 


/  d’^dS-MPA)\A\, 

J  A 


where  |.4|  is  the  “surface  area”,  i.e.,  length,  of  A.  For  Equation  III. 13,  we  use 


[  T]dV  cx  r]{Nv)\Vl 

J  1/ 


where  Nv  is  the  node  associated  with  V  and  |y|  =  fy  dV  is  the  “volume”,  i.e.,  area, 
ofR.[Ref.  1] 

Numerically  evaluating  Equations  III. 12  and  III.13  yields  the  EVE  discretiza¬ 
tion  of  Equation  III. 8.  The  value  associated  with  each  node,  say  (p,  q),  is  a  sum 

ty=p+l,^r=g-f  1 

V  ’I 

W~p—l,z=:q  —  \ 

where  the  coefficients  are  determined  by  the  integration.  A  convenient  way  to 
represent  the  sum  associated  with  each  grid  point  is  to  place  the  coefficients  ayj,z  into 
a  stencil, 

^P+li9 

Ctp^q-l  ap^q  CTp.ij  +  l  '“p,9  ~  ,q  +  <^p,q-\'^p,q-l 


ffi  ^P>q^V>Q  ffi  ,q^p-\-i ,q  ffi  ®P>9+1^P>9+1 
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where,  for  this  example,  the  corner  values  in  the  stencil  are  all  zero.  The  value  for 
node  {p.q)  is  determined  by  applying  the  stencil  to  Up^g. 

To  see  what  type  of  stencils  are  produced,  consider  the  discretization  on  a 
uniform  mesh  with  grid  size  h  =  -^  in  both  coordinates.  We  use  double  subscripts 
p,  (/  =  0  :  m  —  1,  where  0  corresponds  to  the  Neumann  boundary  nodes  and  m  —  1 
corresponds  to  the  Dirichlet  boundary  nodes.  Written  in  stencil  form,  the  interior 
nodes.  0  <  p,  (/  <  m  —  1 ,  for  the  equations  in  Equation  III. 11  are  as  follows  (XI 
indicates  the  sum  of  the  outer  stencil  entries); 


p{ph,iq-l)h)  -E  P{ph,{q  +  l)h) 

p{{p-\)h,qh) 

For  the  corner  Neumann  node  (p,  g  =  0),  the  stencil  is  given  by 


<?/0- 


-E  Mo4) 


h'^ 


h 


h 


“0,0  =  0)  -  -  V’i(0, -)  +  ^^1(7, 0)  . 


For  both  stencils,  the  coefficient  of  the  unknown  to  which  the  stencil  is  applied  is 

-E. 

For  grid  points  adjacent  to  a  Dirichlet  node,  the  stencil  entry  reaching  to  the 
Dirichlet  node  value  is  moved  to  the  right-hand  side  as  a  coefficient  of  the  boundary 
data.  Otherwise,  the  equations  for  these  points  are  the  same  as  above.  For  example, 
at  the  point  {0<p<m  —  l^q  =  m  —  1)  we  have 

p{{p  +  l)h,l  -h) 


p{phAq-l)h)  -E 


p{{p-l)hA  -h) 


Up,m-i  =  h?r]{ph,  I  -  h) 


-p{ph,\  -  -)0o(p/i,l), 


where  E  is  the  sum  of  the  outer  stencil  entries  without  the  boundary  terms  removed, 

=  Piph,  1  -  ■^)  +  {q  -  ^)h)  +  pUp  +  ^)h,  l-h)+  p{{p  -  i)h,  1  -  h). 
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(see  [Ref.  1]) 


2.  Application  to  Axisymmetric  Heat  Equation 

This  technique  is  now  used  to  compute  the  FVE  stencils  for  Equation  II. 2,  by 
combining  the  finite  elements  (see  Figure  2)  with  the  control  volumes  (see  Figure  .3). 
A  control  volume  for  the  Ith.  grid  point,  call  it  V;,  is  a  toroidal  prism  (see  Figure  6). 
It  is  generated  by  taking  the  two-dimensional  sub-volume  for  that  point  (in  the  (r,  ~r) 
|)lane),  and  sweeping  through  360°  in  the  (/?  direction. 

1  z 

^  y 

X  r 


Figure  6.  Toroidal  volume  for  the  conduction  problem. 


Integrating  Equation  III. 4  over  all  control  volumes  Vi,  where  V  —  M 

partitions  the  domain  fl,  we  have 


T^+^dV  =  0  -p  (1  -  TW,  (111.15) 
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or,  upon  application  of  the  Gauss  Divergence  Theorem  the  volume  integral  of  the 
term  becomes  a  surface  integral,  and  we  have 

1  f  -aK  I  -^3  =  -  f  TW  +  (1  -  a)K  [  VT”  •  hdS,  (III.16) 

g  Jv  Js  g  Jv  Js 

where  n  is  the  outward  normal.  Substituting  the  expression  for  approximating  the 
unknown  function  T  (Equation  III. 7),  we  have 


where  we  now  have  integrals  over  each  of  the  (A^  +  1)'^  control  volumes,  resulting  in 
a  set  of  (A^  +  1)^  equations.  This  set  of  equations  can  be  written  both  as  an  operator 
equation,  ==  /(T^),  and  as  a  matrix  equation^,  ~  where  the 

operator  L  and  the  matrix  M  are  given  by 

^ 

M  =  -  /  (t>f+^dV  -aK  [  •  hdS, 

g  -jv  J s 

and  J{T''^)  are  given  by 

/(T")  =  ^  +  (1 TW, 

r  r  1 

fiT^  =  E  -  4>?dV  + {I- a)K  ycf>^-hdS  TP . 

/=i  ^ ^ 

Any  function  whose  coefRcients  satisfy  the  resulting  set  of  discrete  equations  neces¬ 
sarily  satisfies  the  conservation  law  over  any  volume  made  up  of  the  union  of  control 
volumes,  except  possibly  at  the  boundaries.  The  Neumann  conditions  at  r  =  0  and 
2:  =  0  are  incorporated  indirectly  into  T,  by  substituting  the  (zero)  normal  derivative 

^The  operator  may  be  represented  as  a  continuous  operator,  L;  as  a  discrete  operator,  or  as 
a  matrix,  M. 


value  into  the  appropriate  term  in  Equation  III. 15.  The  Dirichlet  conditions  at  the 
I’af  boundaries  are  imposed  directly  on  T\  control  volumes  that  would  usually  be  as¬ 
sociated  with  these  boundary  nodes  are  eliminated  (see  Figure  7).  Thus,  the  reduced 
collection  of  control  volumes  no  longer  partitions  the  grid  17,  which  slightly  impairs 
conservation  in  the  discretization.  However,  since  we  require  the  temperature  to  fall 
off  to  zero  at  these  points,  this  loss  of  conservation  is  not  a  concern.  [Ref.  1] 


0  1  2  N 

r  _ ^ 


Figure  7.  The  problem  domain  cross-section  is  depicted  after  partitioning  into  sub¬ 
volumes.  Homogeneous  Dirichlet  boundary  conditions  obviate  the  necessity  to  define 
volumes  for  grid  points  at  the  far  boundaries. 

The  integrals  over  the  basis  functions  have  been  computed  using  Maple  soft¬ 
ware  and  a  program  designed  by  David  Canright  (see  Appendix  B  and  [Ref.  8]).  In 
order  to  calculate  the  integrals,  we  must  be  able  to  represent  the  unknown  function 
over  the  basis  functions,  which  themselves  are  collections  of  triangular  planes.  There¬ 
fore,  the  method  is  to  determine  first  the  function  for  the  plane  through  three  given 
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|)oint,s  (the  cases  of  vertical  planes  and  three  points  collinear  are  not  considered), 
rids  ])tocedure  is  then  used  to  create  the  triangular  elements  which,  when  assem¬ 
bled,  form  the  hat  function  (see  Figure  5).  For  a  given  grid  point,  say  {k,j),  six  of 
these  triangular  planes  are  joined,  corresponding  to  the  six  triangles  surrounding  the 
point,  NNE,  ENE,  SE,  SSW,WSW,  NW  (see  Figure  4).  The  unknown  function  is 
then  interpolated  over  the  six  elements. 

Once  the  unknown  function  is  represented  by  the  combination  of  these  six 
interpolations,  we  can  use  Maple  to  compute  the  volume  integrals,  and  the  surface 
integrals  of  the  gradients,  in  Equation  III.  17.  The  results  of  these  integrals  provide 
the  coefficients  which  comprise  the  FVE  stencils. 

Thus,  we  have  the  following  stencils  for  the  volume  and  surface  integrals  re¬ 
spectively  for  interior  grid  points: 

32i-5  16; -1-5 

£  =  224;  32;-hll  Tkj,  (III.18) 

p  q  ao4 

16;  -  5  32;  +  5 

and 

f 

i  (EEVT,,A,)<<5  =  |  2j-1  -8j  2j  +  1  nj.  (111.19) 

2; 

where  the  27r  resulting  from  integration  in  the  tp  direction  has  been  factored  out. 
(Recall  that  in  cylindrical  coordinates. 


dV  = 

Iv  Jo 


rdzdrdif.) 

Jr  1  Jz,  1 


Note  that  control  volumes  increase  with  radial  distance  from  the  origin,  which  is 
reflected  in  a  radial  bias  in  the  stencils.  More  specifically,  rj  =  jh  is  the  radial 
distance  from  the  origin,  with  ;'  the  column  index  for  the  unknown  matrix,  h  the 


meshsize.  Thus,  stencil  values  increase  with  distance  from  the  origin. 


These  stencils  are  applied  to  the  unknown  at  a  point  on  the  grid,  designated 
by  its  row/column  position  {k,j).  Stencil  entries  indicate  the  values  to  be  applied; 
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blank  entries  indicate  a  value  of  zero.  Entry  position  indicates  to  which  grid  point 
the  value  is  applied;  left/right  and  up/down  in  the  stencil  correspond  to  neighboring 
pihnts  in  those  directions  on  the  grid.  That  is,  if  stencil  values  are  replaced  by  the 
positions  to  which  they  apply,  we  have 

{k  +  l,j)  (A;  +  1,  j  +  1) 

(kj)  (kj  +  l) 

(k-lj-l)  (k-lj) 

Thus,  for  example,  the  value  that  appears  in  the  [k  +  l,j)  stencil  position  is  applied 
to  the  grid  point  that  occupies  that  same  position. 

Using  techniques  similar  to  those  outlined  in  the  two-dimensional  example, 
boundary  point  stencils  have  been  computed  for  the  volume  and  surface  integrals 
respectively.  i\t  the  origin,  r  =  0,  2  =  0: 

1  5  1 

15  3  To.o,  and  ^  -3  2  To,o  +  j  qir)rdr-  (III.20) 

at  ?•  =  0: 

1  5  1 

h 

25  11  Tfc.o,  and  -  -6  4  (III.21) 

O 

6  [  1 

and  a.t  0  =  0: 

32i-5  16i  +  5.  4i 

24j-8  112i  +  5  8i  +  3  Toj,  and  -  2j  - 1  -8j  2j  + 1  Toj. 

(III.22) 

For  points  adjacent  to  the  far  boundaries,  the  stencils  in  Equations  111.18,  III.19,  III. 21, 
and  1 1 1.22  are  applied.  However,  since  the  homogeneous  Dirichlet  conditions  dictate 
that  these  boundary  values  are  zero,  the  resulting  contribution  after  stencil  appli¬ 
cation  remains  zero.  Thus,  in  effect,  far-boundary  values  do  not  contribute  to  the 
stencils  for  points  adjacent  to  these  boundaries. 


We  can  now  combine  all  of  the  above  stencils  to  generate  the  operator 


=  /‘(T”),  (III.23) 

where  h  is  the  step  size  on  the  grid.  The  matrix  representation  of  M,  is  a  block 
tridiagonal  {N  +  l)'^x(A^  +  1)^  matrix  of  the  form 

A  A  0  0  •••  0  0 

r  A  A  0  •••  0  0 

0  r  A  A  •••  0  0 

M=  0  0  0  : 

:  !  0  i 

0  .  0  r  A  A 

0  .  0  r  A 

where  A,  A,  and  T  are  {N  l)x(A^  +  l)  generic  banded  matrices  (not  all  identical);  A 
is  tridiagonal,  A  is  upper  bidiagonal,  and  F  is  lower  bidiagonal.  When  M  multiplies 
the  matrix  of  [N  +  1)^  values  of  T,  arranged  lexicographically  as  a  column  vector, 
the  matrices  A,  A,  and  F  produce  the  effect  of  the  stencils  “reaching”  function  val¬ 
ues  respectively  on  the  current  row,  the  row  above,  and  the  row  below.  Numerical 
experiments  using  Matlab  to  construct  the  matrix  M  for  N  =  8,12,16,20,  and  24, 
with  g  =  h  and  cr  =  5  and  1,  indicate  that  it  has  full  rank.  Thus,  we  expect  a  unique 
solution  to  the  linear  algebra  problem  which  is  a  discretization  of  the  heat  equation. 
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IV. 


RELAXATION  SCHEMES 


The  FVE  method,  with  weighted-average  time-stepping,  has  been  used  to 
discretize  the  continuum  problem, 

dT 

^  =  /cV^T,  (IV.l) 

=  f\T^)  ■  (IV.2) 

or,  in  matrix  form, 

on  a  grid  of  meshsize  h  =  ^.  The  input  for  these  equations  is  the  solution  at  the 
current  time  step,  T"  (where  T”  and  T”  are  used  interchangeably),  the  result  of 
solving  the  equations  is  the  solution  at  the  next  time  step,  T'^+^,  As  was  indicated 
in  Chapter  III,  there  is  good  reason  to  believe  that  the  matrix  M  is  of  full  rank  and, 
thus,  expect  a  unique  solution  to  exist  for  the  linear  algebra  problem  that  arises  from 
the  discretization  of  the  continuum  problem.  Solution  by  direct  methods  requires 
factorization  of  M  and,  since  M  is  both  large  and  sparse,  solution  by  direct  methods 
may  I')e  impractical.  We  therefore  consider  iterative  methods  to  solve  the  matrix 
ecpiation  at  each  time  step.  These  methods  generate,  for  each  time  step,  a  sequence 
of  approximate  solutions,  {?(«)}';  the  choice  of  relaxation  scheme  determines  whether 
or  not  the  sequence  converges.  Upon  implementation  of  a  relaxation  scheme 

that  produces  a  converging  sequence  of  approximate  solutions  at  each  time  step,  we 
begin  with  an  initial  temperature  distribution,  T^,  and  the  solutions  are  stepped  in 
time.  The  desired  result  is  a  sequence  of  solutions,  corresponding  to  a  sequence  of 
time  steps,  which  tends  to  steady  state.  This  process  allows  for  the  evaluation  of 

^Here  the  subscript  (s)  indexes  the  sequence  of  approximate  solutions;  elsewhere,  subscripts 
without  parentheses  indicate  grid  position  or  vector  components. 
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various  values  for  the  time  step,  as  well  as  of  various  weightings  used  in  averaging  the 
current  and  subsequent  time  steps. 

A.  ITERATION  MATRICES 

It  is  common  in  constructing  iterative  methods  to  propose  a  splitting  for  the 
matrix  M  =  E  —  F,  where  linear  systems  of  the  form  Ex  =  b  are  “easy”  to  solve  (see 
[Ref.  9]).  The  sequence  of  approximations,  (Tis)},  is  generated  by  Er(5+i)  =  FT(s)  +  f 
and,  therefore,  it  is  natural  to  define  an  iteration  matrix,  P  =  E"T,  so  that 
f^s+\)  =  P^(s)  +  E"'/.  Additionally,  if  f*  is  the  exact  solution  so  that  MT*  =  /, 
then  ET*  =  FT*  +  / and  f*  =  E~^Ff*  +  E-^f.  Hence  the  solution,  f*,  is  a  fixed 
point  of  the  iteration  T(s+i)  =  E“TT(5)  +  E“’/*. 

The  matrix  P  =  E~T  is  also  called  the  error  propagation  matrix,  since  if 

=  Pf(,)  +  E-*/ or  =  Pf(,)  +  B  {B  a.  constant  vector),  and  is  the 

error  at  the  (s  +  l)st  step,  then 

rp  rp* 

^(s+1)  “  ^ 

=  [Pf(,)  +  B]-  [Pf*  +  B]  (since  f*  =  FT*  +  B) 

=  Pf(,)-Pf* 

=  p(f(,)-n. 


and  therefore 

?(•+!)  =  P?w-  (IV.4) 

Using  induction,  it  is  easy  to  show  that  €(*,  =  P’e(o).  By  Equation  IV.4,  Cji,  =  Pe(o) 
and 


6(2)  =  Pe(i) 

=  P[Pe(o)] 

=  P^e(o). 
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By  induction,  assume  that  =  P^^(o)i  in  order  to  show  that  6(^+1)  = 


e(i:+l)  =  Pe(A:) 


=  P[P^'e(o)]  (by  the  induction  hypothesis) 

Thus,  since  e^k+i)  =  we  conclude  e^^s)  =  P^^(o)?  for  5  any  integer.  Addi- 


b  ion  ally, 


which  will  be  useful  later. 


ie.ll  =  ||P'e(o)ll  <  l|P!l‘l|e(o)||, 


(1V.5) 


One  of  the  simplest  relaxation  methods  is  the  Jacobi  method,  also  referred 
to  as  simultaneous  displacement  (for  a  more  detailed  account  of  all  of  these  methods, 
see  [Ref.  3]  or  [Ref.  9]).  It  is  produced  by  solving  the  Ith  equation  of  Equation  IV. 3 
for  the  Ith.  unknown,  T;,  /  =  1:(/V+1)^.  Before  proceeding  to  a  discussion  of  the 
iteration  matrix,  we  present  the  component  form  for  this  iteration  scheme; 


rp{new)  _  rp{new) 

^  hj  =  ^  I 


^224i  +  ^8j 


(IV.6) 


where  Ylv  J2s  respectively,  the  sum  of  volume  stencil  entries  and  surface 
stencil  entries  applied  to  the  current  values  of  the  unknown,  except  at  the  grid 
point  on  which  the  stencil  is  centered: 


Yiv  =  +(32i  -  n)T£'?, 


(32j  -  5)r<:'£  +(16j  + 

+(32j  +  ll)T<"> 


(IV.7) 


+( 163  -  +(323 + 


Es  =  +(-i+23)r£">,  +(i  +  2j)rS', 

,  -HjTl-'t]. 


(1V.8) 
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One  iteration  sweep  consists  of  computing  Equation  IV. 6  for  each  component  of  the 
unknown  vector  T.  Provided  that  the  process  converges,  the  sweeps  continue  until  a 
desired  level  of  convergence  is  reached. 

•Another,  more  succinct,  way  to  present  this  method  is  in  matrix  form.  If 
we  write  the  operator  matrix  as  the  sum  of  a  diagonal  matrix  (D),  a  strictly  upper 
triangular  matrix  (U),  and  a  stictly  lower  triangular  matrix  (L), 

M  =  D  +  U  +  L, 

tlie  matrix  equation  to  be  solved  becomes 

(D  +  U  +  L)[f]=/. 

Now  define  E,/  =  D  and  Fj  =  —  (U  +  L),  so  that  this  may  be  written  as 

D[f]  =  _(U  +  L)[f]  +  /,  or 

Ej[f]  =  Fj[f]  +  f, 

or  as 

f  =  -D-*(U  +  L)[f]  +  D-'/,  or 
f  =  E-/Fj[f]  +  E-/f, 

which  corresponds  to  solving  the  /th  equation  for  T;,  /  =  1  ;  [N  +  1)'^.  If  we  define 
the  .Jacobi  iteration  matrix  as 


Pj  =  or 

=  -D-‘(U  +  L), 

the  matrix  form  of  the  .Jacobi  method  becomes 

finew)  ^  ^  Ej^f.  (IV.9) 
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A  modified  version  of  this  method,  called  the  weighted  Jacobi  method,  is 
det.ermioed  by  introducing  a  weight,  0  <  u;  <  1  (co  =  1  is  the  original  Jacobi  method) 
(see  [fief.  3]).  The  matrix  form  is 

]+»<£;'/,  (iv.io) 

where  I  is  the  indentity  matrix,  =  D,  and  =  (1  —  a;  )D  —  cu^TJ  -|~  L),  and 

=  E-'F^ 

=  D-*[(l -co)D-co(U  +  L)] 

=  (1  -  u;)I  +  coPj 


'j^(new)  _  jp 


is  the  weighted  Jacobi  iteration  matrix.  In  this  method,  the  new  approximation  is  a 
weighted  average  of  an  intermediate  approximation  and  the  old  approximation;  the 
intermediate  approximation  is  the  Jacobi  iterate  of  the  old  approximation. 

The  weighted  Jacobi  method  computes  all  of  the  components  of  the  new  ap¬ 
proximation  before  it  begins  to  use  them  in  the  next  approximation.  This  requires 
storing  both  the  current  and  new  approximations;  a  simple  modification  allows  for 
updating  the  current  approximation  “in  place”,  using  only  the  storage  required  for 


one  approximation.  The  Gauss-Seidel  method  incorporates  the  new  changes  as  soon 
as  they  are  computed  by  overwriting  the  current  approximation  with  the  new  (see 
[Ref.  3]).  More  important,  however,  is  that  (on  model  problems)  the  Gauss-Seidel 
method  generally  converges  about  twice  as  fast  at  the  Jacobi  method  (see  [Ref.  9]). 
In  component  form,  this  iteration  scheme  is 


ji{new)  _  rj\{new) 


i224j  +  2 


'(£s)l 


(IV.ll) 


where  the  arrow  indicates  overwriting,  and 


+(32j  -  ll)Ti“";’ 

+(16J  - 


(32j  - 

+(32j  +  S)TtTj 


+(i6j  +  5)r<;w+. 

+(32j  +  iijrgf, 


(IV.12) 
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a.Tjd 

ts=  +(-i+2j)ri“’  +(i+2i)r«  .  (IV.13) 

^  +2jT<!'”'  ) 

Ill  matrix  form,  using  M  =  D  +  U  +  L,  Eg  =  (D  +  L),  and  Fg  =  — U,  we  have 

(D  +  L)|f|  =  -U|f)+/,  or 
Ec[fl  =  Fo(f|  +  /, 

t)r 

finew)  ^  -(D  +  L)-’U[f(°'‘^’]  +  (D  +  L)-7'',  or 

finew)  ^  E5*FG[f(°'‘^’]  +  E7/. 

This  corresponds  to  solving  the  Zth  equation  for  7),  using  the  new  approximatioirs  for 
components  1,2,...,/—  1.  Now,  define  the  Gauss-Seidel  iteration  matrix, 

Pg  =  E^'Fg 

=  -(D  +  L)-'U, 

SO  that  the  iteration  scheme  in  matrix  form  becomes 

finew)  ^  p^[f(old)^  +  p-lf^  (IV.  14) 

With  the  above  lexicographic  ordering  of  the  (A^  +  1)^  components  of  T  and  the 
components  updated  in  ascending  order,  the  effect  of  a  Gauss-Seidel  sweep  is  to  start 
at  ?•  =  0  and  update  in  the  radial  direction  for  each  vertical  step,  starting  at  z  =  0 
(see  Figure  8). 

As  with  the  weighted  Jacobi  method,  we  can  make  a  modification  to  the  Gauss- 
Seidel  method.  Define  a  parameter  7  G  7?.  (7  =  1  is  the  original  Gauss-Seidel  method) 
and  the  modified  method  is  successive  over-relaxation,  SOR,  (see  [Ref.  9])  which, 
in  matrix  form,  is  given  by 

finew)  ^  p^f(old)  ^E^/,  (IV.  15) 


32 


ooooooooo 


^  O  old  approximation 

■  new  approximation 
Figure  8.  Gauss-Seidel  sweep. 

where,  with  =  D  +  7L  and  =  (1  —  7)D  —  7U,  we  have 

P,  =  (D  +  7L)-M(1-7)D-7U] 

= 

Similar  to  the  weighted  Jacobi  method,  SOR  is  a  weighted  average  of  an  intermediate 
approximation  and  the  old  approximation. 

The  question  of  interest  now  is  whether  or  not  the  sequence  of  approximations, 
{7’(s)},  generated  by  ET(s+i)  =  FT(s)  +  /,  converges.  Therefore,  we  make  use  of  the 
following  theorem: 

Theorem  IV.  1  Suppose  that  f  €  7?."  and  M  =  E  —  F  €  is  nonsingular.  If 

E  is  nonsingular  and  the  spectral  radius  of  E“^F  satisfies  p(E“^F)  <  1,  then  the 
iterates  defined  by  ETp+ij  =  FTp)  +  f,  converge  to  T  =  M“^/  for  any  starting 
vector  T(o). 

The  proof  is  found  in  [Ref.  9],  and  makes  use  of  the  following  lemma: 
Lemma  IV. 1  If  p(E~^F)  <  I,  then  (E“^F)^  — >■  0  as  k  00. 
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II 

00 

N  =  12 

=  16 

IV  =  20 

N  =  24 

Jacobi 

.9259 

.9544 

.9675 

.9749 

.9796 

.9721 

.9826 

.9875 

.9903 

.9921 

.9442 

.9652 

.9750 

.9806 

.9841 

Gauss-Seidel 

.8395 

.8982 

.9263 

.9425 

.9530 

7  =  S 

.9670 

.9793 

.9851 

.9884 

.9905 

7=i _ 

.9189 

.9488 

.9630 

.9711 

.9764 

Table  I.  Spectral  Radii  for  Crank-Nicholson  Time  Stepping  {a  =  |). 


N  =  S 

N  =  12 

N  =  16 

=  20 

II 

Jacobi 

.9499 

.9709 

.9801 

.9850 

.9881 

-=3  . 

.9817 

.9892 

.9925 

.9943 

.9955 

.9634 

.9784 

.9850 

.9887 

.9909 

Gauss-Seidel 

.8931 

.9362 

.9556 

.9663 

.9730 

7  =  ^ 

.9782 

.9871 

.9910 

.9932 

.9946 

7  =  ^ 

.9462 

.9680 

.9777 

.9831 

.9865 

Table  II.  Spectral  Radii  for  Implicit  Time  Stepping  (cv  =  1). 

The  matrices  M,  Ej,  and  Eg,  and  iteration  matrices,  P,  for  various  values 
of  N ,  a  =  ^  and  1,  and  g  =  h,  have  been  constructed  using  Matlab.  We  have 
experimentally  verified  that  M,  Ey,  and  Eg  are  nonsingular,  and  the  spectral  radii  of 
the  error  propagation  matrices,  P,  have  been  computed.  Due  to  memory  limitations 
with  Matlab,  the  calculations  are  made  for  relatively  small  values  of  N.  The  results 
are  indicated  in  Tables  I  and  II,  where  the  value  of  a  indicates  the  type  of  time 
stepping,  and  the  values  of  u  and  7  are  the  weightings  in  the  weighted  .lacobi  and 
SOR  methods  respectively. 

The  results  of  these  numerical  experiments  indicate  that  both  /)(Pj)  <  1  and 
p(Pg)  <  1,  with  p(Pg)  <  p(Pj)-  The  modifications  to  the  Jacobi  and  Gauss-Seidel 
methods  do  not  provide  any  improvement;  the  spectral  radii  for  both  methods  are 
higher  for  the  modified  iteration  matrices  than  for  the  original  iteration  matrices. 
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Thus,  the  Gauss-Seidel  method,  used  with  either  the  Crank-Nicholson  or  implicit 
time-stepping  scheme,  appears  to  be  the  best  choice  of  these  relaxation  schemes.  How¬ 
ever,  if  the  spectral  radius  of  P  is  near  unity,  convergence  may  by  unacceptably  slow 
since  the  error  tends  to  0  like  /9(P)^,  as  indicated  by  the  lemma  and  Pll  <  mi'lirmll 
(Equation  IV. 5).  As  is  evident  from  the  above  tables,  even  for  a  moderately  spaced 
grid  (e.g.,  N  =  16,  g  —  h,  and  a  =  |  or  1),  p(P)  >  .9,  and  the  spectral  radius 
increases  with  N. 

B.  ALTERNATIVE  SCHEMES:  LINE  RELAXATION 

Both  the  .Jacobi  and  Gauss-Seidel  methods  give  rise  to  iteration  matrices,  are 
implemented  in  a  straightforward  manner,  and  are  attractive  because  of  their  sim¬ 
plicity.  While  the  ease  of  implementation  is  a  significant  advantage,  the  convergence 
properties  of  either  or  both  may  not  be  acceptable  and,  therefore,  alternative  schemes 
may  be  desirable.  One  such,  which  seems  reasonable  based  on  the  geometry  of  the 
lu'oblem,  is  line  relaxation.  There  are  two  obvious  options  in  this  regard:  radial  or 
vertical  line  relaxation  (see  Figures  9  and  10),  where  either  an  entire  row  or  column 
is  updated  simultaneously.  Both  options  require  the  solution  of  an  {N  -f  l)x(iV  +  1) 
tridiagonal  matrix  for  each  row/column  update  in  the  unknown  matrix.  That  is,  while 
the  previously  outlined  relaxation  schemes  (point  relaxation)  proceed  by  successively 
solving  a  collection  of  algebraic  equations  for  one  unknown,  line  relaxation  proceeds 
fry  solving  a  succession  of  matrix  equations.  For  example,  suppose  that  the  matrices 
.4,  A,  and  F  are  tridiagonal,  upper  bidiagonal,  and  lower  bidiagonal,  respectively,  and 
that  Ti  and  fi  are  the  portions  of  the  respective  unknown  and  right  hand  side  vectors 
in  Equation  IV. 3  that  are  rows  in  their  corresponding  matrices: 
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R  adial  line  relaxation  consists  of  solving  the  following  succession  of  matrix  equations: 


Ti  f-  A~\f,-AT2), 

T2  A-\f2-TT,-AT:,l 
Ts  e-  .4->(/3-rr,- Ar4),  and 
T,  e-  A-\U-rT3). 


The  computational  cost  of  updating  an  entire  row/column  at  a  time  is  higher 
than  a  row/column  update  using  either  the  Jacobi  or  Gauss-Seidel  methods.  However, 
this  type  of  relaxation  may  be  sufficiently  efficacious  to  warrant  the  extra  expense,  in 
that  convergence  may  be  achieved  significantly  faster.  For  comparison,  we  now  con¬ 
sider  local  mode  analysis  of  the  Jacobi  and  Gauss-Seidel  methods  and  line  relaxation. 


latest  column  update 


next  column  update 


I  i 


r 


“  "  —  old  approximation 
new  approximation 

Figure  10.  Updating  an  entire  column  at  a  time. 

C.  LOCAL  MODE  ANALYSIS 

While  an  examination  of  the  spectral  radii  of  iteration  (error  progagation) 
matrices  can  be  instructive,  it  is  often  useful  to  conduct  a  more  detailed  analysis. 
What  follows  is  a  local  mode  analysis  of  various  iteration  schemes.  Since  the  error 
propagation  matrices  indicate  how  the  error  evolves  during  the  iteration  process,  we 
use  a  DFT  (see  Appendix  A,  [Ref.  10])  to  expand  components  of  the  error  equations 
associated  with  the  various  relaxation  schemes.  The  coefficients  in  these  expansions 
are  the  factors  by  which  the  corresponding  modes  of  the  error  are  magnified/reduced 
for  each  relaxation  sweep.  Thus,  by  examining  these  coefficients,  we  can  determine 
how  quickly  the  error  tends  to  zero,  which  indicates  how  quickly  the  solution  con¬ 
verges.  The  following  analysis  does  not  apply  to  points  on  the  boundary,  it  is  only 
valid  for  interior  points.  Thus,  while  it  gives  more  information  about  the  functioning 
of  the  scheme  in  the  interior  of  the  domain,  boundary  peculiarities  are  not  addressed. 


38 


We  begin  by  recalling  that  the  current  error,  6(5),  is  the  difference  between 
i.lie  exact  solution,  T"',  and  the  current  approximation, 

=  r*  —  Tr.ou 


where  we  desire  that  the  sequence,  {^(5)},  tend  to  zero.  Substituting  this  expression 
into  Equations  IV. 6  and  IV.  11,  we  have  the  following  error  equations  for  the  Jacobi 
and  Gauss-Seidel  methods  where,  in  order  to  avoid  confusion  with  the  exponential 
Function  in  the  DFT  expansion,  we  represent  the  components  of  e  by 


Jacobi 


Seidel 


(new) 

'■ij 


)  j-12^5  J 

J-224i  +  ‘f!8i 

^224i  +  ^8j  ’ 


(IV.22) 


(IV.23) 


where  £5  are  now  applied  to  the  (see  Equations  IV. 7,  IV. 8,  IV.  12, 

and  IV,  13).  Equations  IV.22  and  IV.23  indicate  how  the  sequence  {6(5)}  evolves 
during  the  iteration  process.  If  we  expand  these  relationships  in  a  discrete  Fourier 
transform,  (DFT)  (see  Appendix  A  or  [Ref.  10]),  the  magnitude  of  the  transform  co¬ 
efficients  will  indicate  the  “amount”  of  each  mode  of  the  error  that  is  present  in  each 
component  of  6(5).  We  then  compare  coefficients  to  determine  the  ratio  between  the 
new  and  current  approximations  for  each  component  of  the  error.  In  other  words,  the 
DFT  allows  us  to  analyze  the  growth  of  the  error  by  examining  component  behavior. 

Expanding  Equation  IV.22  in  a  DFT,  where 


m=-f+l 


ilnkl  i2irjm 

Vi  me  N  e  N  _ 


C{l,m) 


t2nkl  i2njm 

e  N  e  N  0,  and  J~\3)  = 


^224i  +  ^8i 


and  the  superscripts  (o)  and  (n)  indicate  the  components  of  the  old  and  new  approx¬ 
imations,  we  have 


+(16,+5)v£>C(/,m)e^ 


+(32i  -  ll)K^C(/,77i)e 


+(32j  +  ll)<Ja(/,m)e"^ 


+(16i-5)V,;,;’C((,m)e - ^ 

+(32i  +  5)V^2C(/,m)e-‘^) 

O^l^h  ,  ,  iZrrl 

- ^(2jV;JJC(/,m)e  N 

+{-1  +  2j)Vi[°lCil,  m)e~'^ 
+(1  +  2j)v!'°lC{l,in)e'^ 
+2iV;2C(/,m)e^^)  ]  }T{j), 


Making  use  of  the  orthogonality  property  of  the  complex  exponential  (see  Appendix  A) 
we  can  equate  individual  terms  of  the  sum,  and  then  divide  by  C{l,m)  to  give 

vS  =  -  j^((32j-5)V;;,:'eT+(16i  +  5)ViSe'^ 
+(32i-ll)V;,l:'e=^+(32j  +  llWj:'e'^ 
+(16i-5)l^I:'e-'^+(32i  +  5)V5;:'e-'^) 


+(H-2i)V^,<:>e‘^+2iVi;:V 


‘)]  m, 


.  .  ■2]ri  .  .  .  t27r(l+: 

— — ((32;  -  5)e  TV  -p(i6;+5)e  tv 


+(32j-ll)e^ 


^  ^  \  *27rm 

+  (32;  +  ll)e  tv 


t27r(f-hm)  t2nl  . 

+(16;  —  5)e  tv  q.  (32_^  +  5)e  tv  ) 

OCKh  .  t27rf  /  /-V  •\  t27rm 
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N  =  16 

N  =  32 

J  =  1 

■  _  iV  ■' 

J  2 

j  =  N-l 

J  =  1 

7^ 

j  =  N-l 

a  —  Ty 

1.009 

1.009 

1.009 

1.026 

1.026 

1.026 

a  =  1 

1.005 

1.005 

1.005 

1.013 

1.013 

1.013 

Teible  III.  Maximum  Values  of 
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for  /,m  =  0  :  N  for  the  Jacobi  Method  of 
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Merc 


fore,  the  ratio  of  the  new  coefficient  to  the  old  is 
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In  order  to  determine  the  greatest  factor  by  which  a  mode  of  the  error  is  multiplied, 


we 


seek  the  maximum  of 


1^ 


over  the  values  /,  m=:— ^  +  1:^.  In  other  words,  we 


seek  to  determine  a  bound  on  how  well  we  can  expect  this  type  of  relaxation  scheme 
to  perform.  This  ratio  is  a  function  of  N,  /,  rn  and  j  and,  while  difficult  to  determine 
analytically,  may  be  calculated  numerically.  Using  Matlab,  the  maximum  of  this  ratio 
for  l,7n  =  0  :  N  has  been  calculated  for  N  =  16,32;  a  =  i,  1;  and  j  =  1,  y,  N  —  1;  the 
results  are  indicated  in  Table  III  (see  Appendix  A  for  a  discussion  of  the  equivalence 
of  centered  indices,  /,m  =  —  y  +  1  :  y,  and  non-centered  indices,  l,m  —  0  :  N). 
Additionally,  by  considering  the  matrix  of  grid  points  /,  m  =  0  :  A^,  we  can  determine 
a  correspondence  between  sample  points  on  the  grid  and  type  of  associated  frequency 

|V‘(”)  I 

(see  Figure  11  and  Appendix  A).  The  matrices  of  values  of  for  j  =  N  —  I  are 

IU,ml 

depicted  in  Figures  12  and  13. 


The  information  in  Table  III  indicates  that  the  Jacobi  method  results  in  a 
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maximum  of  greater  than  one;  this  value  does  not  appear  to  depend  on  grid 

position.  In  other  words,  there  is  at  least  some  component  of  the  error  that  is  magni- 
hed.  instead  of  reduced,  by  this  iteration  scheme.  Moreover,  it  appears  that  the  value 
increases  with  the  number  of  intervals  on  the  grid.  Additionally,  Figures  12  and  13 
indicate  that  the  maximum  of  this  ratio  occurs  for  both  low  and  high  frequencies. 
Thus,  it  appears  that  the  .Jacobi  method  will  not  be  effective  in  generating  a  solution 
lo  l  lie  point-source  problem. 


Low  frequency 

Mixed 

Frequency 

Low  frequency 

N 

_  Mixed 

Frequently 

Mixed 

2 

Frequency 

Frequency 

Mixed 

Low  frequency 

Frequency 

Low  frequency 

0 

- i - : 

Figure  11.  A  two-dimensional  representation  of  frequencies  corresponding  to  a  single 
set  of  non-centered  sample  points. 


In  order  to  compute  the  ratio 


for  the  Gauss-Seidel  relaxation  scheme,  we 


rewrite  Equation  IV. 23  as 


(new)  _ 


-  cvi-jr! 

-  Mit  -  GoKll+i) 


(IV.25) 
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1.2-1 
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Figure.  12.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Jacobi  relaxation: 
N  =  32,  a  =  ^,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 


where  A  =  16j-5,  B  =  32i-h5,  C  =  32i-ll,  D  =  224i,  E  =  32i-Ml,  F  =  32j-5,  G  = 
16j/'  -h  5,  and  H  =  2j,  I  =  —  1  -f  2j,  J  =  —8j,  K  =  I  +  2j,  L  =  2j.  Expanding  in  a 
DFT,  equating  terms  in  the  (double)  sum  by  virtue  of  the  orthogonality  property  of 
the  complex  exponential  (see  Appendix  A),  and  dividing  by  we  have 
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Figure  13.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Jacobi  relaxation: 
;V  =  32,  a-  =  1,  j  =  —  1.  Frequency  ranges  are  as  in  Figure  11. 
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Thus,  the  ratio  is  given  by 
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T'  U  AT  A/ 

.As  above,  we  now  seek  to  maximize  over  the  values  /,m  =  — -j  +  1  :  y  in  order 
to  determine  a  bound  on  how  well  this  scheme  performs.  We  again  use  Matlab  to 
compute  the  maximum  of  this  ratio  for  l,m  =  0  :  N  for  N  =  16,32;  a  =  1;  and 

j  —  1,  y,  —  1;  the  results  are  indicated  in  Table  IV.  Additionally,  the  matrices  of 

lyC”)! 

values  of  for  j  =  —  1  are  depicted  in  Figures  14  and  15. 

The  information  in  Table  IV  indicates  that  the  Gauss-Seidel  method  results 
in  a  maximum  of  less  than  one.  Thus,  all  frequency  components  of  the  error 

are  reduced  by  this  relaxation  scheme,  which  is  what  we  seek.  However,  the  value 
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N  =  16 

N  =  32 

i  =  T 

j  =  N-  1 

J  =  l 

i  =  T 

j  =  N-l 

a;  =  ^ 

.8986 

.  8796 

.8781 

.9152 

.8977 

.8970 

a  =  1 

_ 1 

.9478 

.9374 

.9366 

.9567 

.9473 

.9469 

I 


Table  IV.  Maximum  Values  of  7;^  for  l,rn  —  0  :  N  for  the  Gauss-Seidel  Method  of 


Relaxation. 


increases  with  grid  size  N,  and  also  as  the  time  stepping  is  shifted  from  Crank- 
.Ticliolson  to  fully  implicit.  Additionally,  while  all  of  the  values  are  less  than  one, 
as  we  move  to  grids  with  larger  N  and  toward  a  fully  implicit  time  scheme,  they 
l)ecome  close  to  unity.  This  means  that,  although  we  may  reasonably  expect  to 
see  this  relaxation  process  converge,  it  may  be  quite  slow.  Another  point  is  that 
the  maximum  value  depends  on  grid  position;  the  shorter  the  radial  distance,  the 
larger  the  maximum,  which  apparently  is  a  reflection  of  the  radial  bias  of  the  stencils. 
.Additionally,  Figures  14  and  15  indicate  that  the  maximum  of  this  ratio  occurs  only 
over  the  low  frequencies,  and  that  the  values  associated  with  error  components  over 
the  high  frequencies  are  quite  low.  This  performance  over  the  high  frequencies  will 
be  of  importance  in  Chapter  VI.  Thus,  it  appears  that  the  Gauss-Seidel  method  will 
be  effective  in  generating  an  iterative  solution  to  the  point-source  problem,  albeit 
potentially  slow  to  converge. 

We  now  apply  similar  techniques  to  the  radial/vertical  line  relaxation  schemes. 
The  error  equations  come  from  rewriting  Equations  IV. 20  and  IV. 21  respectively  as 
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Figure  14.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Gauss-Seidel  relax¬ 
ation:  N  =  32,  a  "  j  =  N  --  1.  Frequency  ranges  are  as  in  Figure  11. 
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Figure  15.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Gauss-Seidel  relax¬ 
ation:  N  =  32,  a  =  I,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 
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N  =  16 

N  =  32 

3  =  1 

i  =  T 

3  =  N-\ 

J  =  1 

3  =  F 

j  =  N-\ 

.7689 

.7723 

.7726 

.8011 

.8051 

.8052 

a  =  1 

.8758 

.8768 

.8768 

.8954 

.8965 

.8965 

Table  V.  Maximum  Values  of  '/’a  for  l.rn  =  0  :  N  for  Radial  Line  Relaxation. 

"bus,  the  numerators  of  the  ratio  are  given  by 


/i"  .  i2-7rt  i2n-(i+m)  .  OCKiH  .  i2wl , 

-Fe  N  -  Ge  ^ - \-Iq  n  ] 


and 


5r384 


r  i2irTn  i27r(i-f  m)  OLKjH  ^  t2nm  - 

r/384^“^^  ^  -G'e  ^  ]-— [-/Te  ^  ], 


and  the  denominators  are  given  by  . 
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As  above,  we  now  seek  to  maximize 


i<’i 


over  the  values  l,m  =  — ^  +  1  ;  We 


\ya\  ' 

continue  to  use  Matlab  to  compute  the  maximum  of  these  ratios  for  /,m  =  0  :  AT  for 

i'V  =  16,32;  q:  =  0,  1;  and  j  =  1,  y,  —  1;  the  results  are  indicated  in  Tables  V 

and  VI.  Additionally,  the  matrices  of  values  of  for  j  =  A^  —  1  are  depicted  in 

U/,ml 

Figures  16,  17,  18,  and  19. 

The  information  in  Tables  V  and  VI  indicates  that  both  line  relaxation  meth- 
ods  result  in  a  maximum  of  less  than  one;  radial  line  relaxation  seems  to  promise 
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N  =  16 

N  =  32 

j  =  1 

7^ 

j  =  N-l 

j  =  1 

J  =  Y 

j  =  N-l 

.8379 

.7837 

.7787 

.8627 

.8103 

.8079 

a  =  1 

.9149 

.8834 

.8805 

.9290 

.8994 

.8981 

Table  VI.  Maximum  Values  of  — ^  for  /,m  =  0  :  for  Vertical  Line  Relaxation. 

l)etter  performance.  Many  of  the  results  for  the  line  relaxation  schemes  parallel  those 
noted  for  the  Gauss-Seidel  method.  All  frequency  components  of  the  error  are  reduced 
l)y  the  line  relaxation  schemes;  their  performance  would  appear  to  be  slightly  better 
than  that  for  the  Gauss-Seidel  method.  Additionally,  the  maximum  values  increase 
with  both  the  grid  size  N,  and  as  the  time  stepping  is  shifted  from  Crank-Nicholson 
to  fidly  implicit.  As  with  the  Gauss-Seidel  method,  while  all  of  the  values  are  less 
than  one,  as  we  move  to  grids  with  larger  N  and  toward  a  fully  implicit  time  scheme, 
they  become  close  to  unity.  This  means  that,  although  we  may  reasonably  expect 
to  see  this  relaxation  process  converge,  it  may  be  quite  slow,  and  that  the  increased 
amount  of  computational  work  as  compared  with  work  for  the  Gauss-Seidel  method 
may  not  be  worth  the  payoff.  The  maximum  values  for  these  schemes  also  depend  on 
grid  position.  For  the  horizontal  line  relaxation,  the  shorter  the  radial  distance,  the 
smaller  the  maximum;  for  the  vertical  line  relaxation,  the  shorter  the  radial  distance, 
the  larger  the  maximum.  Once  again,  apparently,  this  is  a  reflection  of  the  radial  bias 
of  the  stencils.  Additionally,  Figures  16,  17,  18,  and  19  indicate  that  the  maximum 
of  these  ratios  occurs  only  over  low  frequencies,  and  that  the  values  associated  with 
error  components  over  the  high  frequencies  are  quite  low.  In  fact,  the  values  for  the 
line  relaxation  schemes  over  the  high  frequencies  appear  to  be  about  one-half  the  cor¬ 
responding  values  for  the  Gauss-Seidel  method.  As  noted  earlier,  this  performance 
over  the  high  frequencies  will  be  of  importance  in  Chapter  VI. 

Some  general  trends  are  evident  from  the  results  of  the  numerical  experiments. 
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Figure  16.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  radial  line  relaxation: 
/V  =  32,  a  =  j  =  N  —  Frequency  ranges  are  as  in  Figure  11. 

.Since  the  maxima  should  all  be  less  than  one  to  guarantee  convergence,  it  appears  that 
the  .Jacobi  method  will  not  be  useful  in  iteratively  solving  the  point-source  problem 
(Equation  11.2).  Of  the  schemes  that  have  maxima  less  than  one,  those  that  have  the 
smallest  maxima  are  those  that  reduce  error  components  most  quickly  and,  therefore, 
are  the  schemes  that  will  converge  most  quickly.  The  line  relaxation  schemes  have  the 
smallest  maxima;  the  maximum  values  over  the  high  frequencies  are  about  one-half 
the  corresponding  values  for  the  Gauss-Seidel  method,  but  there  is  little  increase  in 
l)erformance  over  the  low  frequencies.  To  determine  whether  or  the  extra  work  in 
using  line  relaxation  is  worth  the  effort,  a  cost  comparison  of  getting  to  the  same 
level  of  error  as  other  relaxation  schemes  will  need  to  be  done.  Additionally,  the 
maximum  values  vary  depending  on  the  time  stepping  scheme  and  the  grid  position; 
they  are  higher  for  implicit  time-stepping  than  for  Crank-Nicholson,  indicating  that 
convergence  will  take  longer  for  the  former  as  compared  to  the  latter,  and  the  maxi¬ 
mum  values  increase  with  N.  Hence,  as  we  move  to  a  grid  with  more  nodes,  we  not 
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Figure  17.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  radial  line  relaxation: 
/V  =  32,  a  =  \,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 

only  have  more  work  to  do  by  virtue  of  the  greater  number  of  equations  to  solve,  we 
also  have  to  work  harder  to  reduce  each  component  of  the  error.  That  is,  the  effect 
of  moving  to  a  finer  grid  is  to  more  than  quadruple  the  work  required  to  solve  the 
problem. 

The  Gauss-Seidel  and  line  relaxation  schemes  appear  to  be  the  most  efficient 
in  terms  of  reducing  error  components.  Additionally,  in  Chapter  VI,  the  performance 
of  these  schemes  over  the  high  frequencies  will  be  of  interest;  we  will  want  to  use  those 
schemes  that  have  the  smallest  maximum  value  over  the  high  frequencies.  In  particu¬ 
lar,  for  the  Gauss-Seidel  and  line  relaxation  schemes,  the  high  frequency  components 
of  the  error  are  eliminated  much  more  quickly  than  the  low  frequency  components. 
Considering  the  amount  of  work  to  be  done  and  the  performance  of  the  scheme,  we 
implement  the  Gauss-Seidel  method  in  the  solution  process. 
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Figure  18.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  vertical  line  relax¬ 
ation:  N  —  Z2,  a  =  j  =  N  —  Frequency  ranges  are  as  in  Figure  11. 


Figure  19.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  vertical  line  relax¬ 
ation:  N  =  Z2,  a  =  \,  j  —  N  —  Frequency  ranges  are  as  in  Figure  11. 
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V. 


ITERATIVE  SOLUTION 


The  results  of  Chapter  IV  indicate  that  the  Gauss-Seidel  method  is  the  sim¬ 
plest  of  the  relaxation  schemes  that  appear  to  result  in  a  converging  sequence  of 
approximate  solutions,  and  therefore  Gauss-Seidel  is  our  method  of  choice.  Having 
chosen  a  relaxation  scheme,  we  must  now,  as  part  of  the  iterative  solution  process, 
make  suitable  choices  for  the  type  of  time-stepping,  a,  and  the  time  step  size,  g;  the 
goal  is  to  maximize  accuracy  and  minimize  computational  work.  We  now  conduct 
experiments  to  solve  Equation  IV. 3  by  iterating  with  a  solver  which  is  based  on  the 
iteration  matrix.  Pc;  =  — (D  -|-  L)“^U,  and  use  this  process  to  evaluate  some  of  the 
l>ossible  choices  for  a  and  g.  The  technique  is  to  solve  first  at  a  single  time  step  by 
iterating  until  the  difference  between  successive  approximations,  —  T^l)\\  or, 

perhaps  better,  —  /(T(”))||,  is  negligible.  This  approximate  solution,  T",  is 

used  as  input  for  the  time-stepping  scheme  so  that 

1  r  f 

f{T^)  =J2[-  +  (1  -  •  f>'dS]Tr 

jri  g  Jv  Js 

(as  in  Chapter  III)  becomes  the  new  right  hand  side  in  Equation  IV. 3.  The  iteration 
process  is  repeated  to  obtain  the  approximation  at  the  new  time  step.  The  solution 
is  stepped  in  time  in  this  fashion  until  the  difference  between  solutions  at  successive 
time  steps  is  negligible,  representing  a  steady  solution. 

An  algorithm  for  this  process  might  look  something  like  the  following: 
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Given  an  initial  guess,  T°,  and  tolerances  (^i,  (^2, 
Set  n  —  0 

While  ||f"+'  -  f"||  >  5-2 

1)  Compute  /(T”) 

2)  Set  .S’  =  1,  =  f" 

3)  While  |j7;”|'  -  >  3, 

r,’21)  ^  +  (d  +  L)-'/(f”) 

.s  =  s  +  1 
End  while 


4)  ^  f’;+‘ 

Return  as  the  solution  at  t  =  (n  +  [)At 

n  =  n  +  I 
End  while 


^steady  ^  ^ 


n 


Stop 


Starting  with  an  initial  temperature  distribution  of  T  =  0  everywhere,  bound¬ 
ary  conditions  specified  by  Equations  11.18,  11.19,  and  11.20,  N  =  16,  and  a  =  |, 
several  values  of  g  are  used  in  an  attempt  to  achieve  a  reasonable  solution.  The  re¬ 
sults  of  these  experiments  indicate  that  for  larger  values  of  g,  say  g  =  h,  the  solution 
exhibits  a  non-physical  oscillation  in  time  (see  Figure  20)  at  the  origin,  where  the 
heat  source  is  located.  Instead,  the  solution  should  monotonically  increase  until  it 
levels  off  at  steady  state. 

The  oscillation  can  be  removed  by  making  the  temporal  step  size  smaller,  say 
g  =  (see  Figure  21).  This  step  size,  however,  results  in  negative  temperatures 
along  the  z  =  0  axis  near  the  origin.  Negative  temperatures,  like  the  oscillations,  are 
physically  meaningless.  A  valid  solution  process  must  eliminate  both  of  these  non¬ 
physical  characteristics  from  the  solution.  However,  it  turns  out  that  for  cv  =  ^,  there 
is  no  value  of  g  that  produces  a  physically  meaningful  solution.  In  particular,  for 
g  =  /),§ ,  there  are  both  oscillations  and  negative  values  in  the  approximate  solution 
(see  Figure  22  for  a  comparison  of  several  values  of  g).  Thus,  in  order  to  compute  a 
realistic  solution,  it  is  necessary  to  use  values  of  a  >  |. 
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Maximum  Solution  Values 


Figure  20.  The  maximum  values  of  the  solution  for  each  time  step,  where  «  =  i  and 
g  =  h. 


In  an  attempt  to  produce  a  solution  that  is  physically  meaningful,  numerical 
experiments  are  conducted  using  values  of  |  <  o;  <  1,  and  values  of  g  ranging  from 
h?  to  h  for  each  value  of  a.  Solutions  with  no  oscillation  and  no  negative  values  are 
possible  for,  say  a  =  |,  but  this  requires  that  the  temporal  step  size  be  as  small  as 
g  =  The  size  of  the  temporal  step  indicates  how  much  work  will  have  to  be 

done  to  reach  a  steady  solution;  the  smaller  the  step,  the  longer  to  reach  steady  state. 
In  order  to  minimize  the  work,  we  must  maximize  the  temporal  step  size.  Further 
experimentation  indicates  that  for  a  =  1  (fully  implicit  time  stepping)  and  g  =  h, 
the  result  is  not  only  a  physically  meaningful  solution  (i.e.,  no  oscillations  and  no 
negative  values),  but  also  a  reasonable  temporal  step  size.  Therefore,  the  iterative 
solution  is  computed  using  these  values. 

Iterative  solutions  for  N  =  8, 16,  and  32,  are  computed  using  the  same  initial 
temperature  distribution  and  boundary  conditions  as  above.  One  measure  of  accu¬ 
racy  of  these  solutions  is  determined  by  comparing  them  with  an  analytic  solution 
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Maximum  Solution  Values 


Figure  21.  The  maximum  values  of  the  solution  for  each  time  step,  where  cv  =  |  and 
<-/  =  . 

(Equations  11.15  and  11.16).  However,  the  problem  that  gives  rise  to  the  analytic  so¬ 
lution  presented  in  Chapter  II  is  somewhat  different  from  the  linear  algebra  problem 
that  we  are  solving.  Therefore,  as  in  [Ref.  11]  we  approximate  the  exact  solution  of 
Equation  II. 2  by  solving  Equation  IV. 3  for  V  =  64  (see  Figure  23).  We  then  compare 
this  solution  with  solutions  of  Equation  IV. 3  on  coarser  grids  {N  =  32  and  N  =  16) 
to  get  approximations  of  the  discretization  errors,  which  may  be  used  to  give  an  in¬ 
dication  of  the  order  of  the  accuracy  of  the  solution.  The  measure  that  we  use  is  the 
discrete  energy  norm,  defined  by 

1110*111  =  (V.l) 

where  (•,  •)  denotes  the  Euclidean  inner  product,  is  the  operator  defined  in  Equa¬ 
tion  IV. 2,  and  =  T’’'  —  T*  approximates  the  discretization  error,  where  is  the 
solution  to  Equation  IV. 2  on  grid  /i,  and  T*  is  the  “exact”  solution  (i.e.,  solution  of 
Equation  IV. 2  with  iV  =  64  sampled  on  grid  h).  (see  [Ref.  11]) 
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90 


Maximum  Solution  Values 


Figure  22.  The  maximum  values  of  the  solution  for  each  time  step,  where  a  =  “o” 

indicates  g  =  indicates  g  =  h},  indicates  g  =  hi ^  indicates  g  = 

(  ’are  must  be  taken  in  analyzing  the  figure  since  the  apparent  equivalence  of  time 
steps  can  be  misleading,  e.g.,  for  g  =  h^  and  h  =  -T,  it  takes  16  time  steps  to  ‘‘equal” 
one  time  step  for  g  =  h. 


The  solutions  for  an  inital  time  step\  for  a  time  of  one  second,  and  for  steady 
state  have  been  compared  using  the  discrete  energy  norm,  with  the  resulting  dis¬ 
cretization  errors  indicated  in  Table  VII.  We  might  hope  to  see  a  common  factor  of 
decrease  as  we  move  to  finer  grids.  That  is,  if  the  discretization  error  on  grid  h  were 
of  the  order  e  =  for  some  constant  c,  then  the  ratio  of  errors  for  grids  h  and  | 


^  With  g  =  h,  the  size  of  a  time  step  differs  with  grid  size.  We  have  adjusted  for  this  difference  by 
requiring  that  the  same  ‘‘amount”  of  time  elapse  for  solutions  used  in  the  comparison,  e.g.,  for  the 
initial  time  step,  the  “exact”  solution  on  iV  =  64  after  eight  time  steps  is  compared  to  the  solution 
on  N  =  32  after  four  time  steps,  to  the  solution  on  N  =  IQ  after  two  time  steps,  and  to  the  solution 
on  yV  =  8  after  one  time  step. 
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Steady  Solution  for  N=64 


would  be 

chP 

2-p|. 

Thus,  a  decrease  of  a  factor  of  four  would  indicate  that  the  process  achieves  0{h^) 
discretization  error  (see  [Ref.  11]).  Table  VII  indicates  that  the  discretization  error 
decreases  by  a  factor  of  about  two  when  moving  from  grid  size  N  =  8  to  N  =  16, 
and  by  a  factor  of  about  four  when  moving  from  grid  N  =  16  to  N  =  32.  A  grid 
size  N  =  8  may  be  too  small  to  adequately  reflect  the  accuracy  of  the  solution,  which 
would  leave  the  factor  of  four,  possibly  suggesting  that  the  discretization  error  is 
0{h‘^).  Even  so,  a  specific  estimate  of  the  discretization  error  is  best  deferred  until 
more  information  can  be  obtained,  i.e.,  solutions  on  finer  grids. 


N  =  S 

N  =  16 

N  =  32 

Initial  time  step. 

43.07 

24.02 

6.85 

Time  =  “1  second” 

6.87 

Steady  state 

42.99 

23.98 

6.87 

Table  VII.  Discretization  Errors  for  the  Iterative  Solution. 


VI. 


MULTILEVEL  SOLUTION 


.^Jow  that  we  can  solve  a  discrete  representation  of  the  heat  equation  (Equa¬ 
tion  IV. 2)  iteratively,  we  want  to  apply  multigrid  techniques  in  order  to  attempt  to 
accelerate  the  solution  process.  Multigrid  is  a  method  to  improve  on  solution  by 
relaxation  by  making  use  of  the  advantages  of  working  on  successively  coarser  grids 
(for  a  more  complete  treatment,  see  [Ref.  3]).  It  has  been  used  with  marked  suc¬ 
cess  in  solving  a  variety  of  problems,  specifically  elliptic  partial  differential  equations. 
.Altliough  the  problem  we  are  solving  is  parabolic,  in  the  time-stepping  regime  we 
must  solve  an  elliptic  problem  at  each  time  step.  Thus,  multigrid  may  prove  useful 
in  streamlining  our  solution  process.  We  begin  by  introducing  the  multigrid  method 
and  some  of  the  exigencies  of  the  use  of  multiple  grids.  We  then  analyze  some  of  the 
characteristics  of  the  method  using  local  mode  analysis,  and  present  numerical  re¬ 
sults  from  the  solution  process.  The  amount  of  work  to  achieve  the  iterative  solution 
serves  as  a  baseline  against  which  the  computational  work  for  the  multilevel  solution 
is  compared. 

As  indicated  in  the  analysis  of  the  Gauss-Seidel  and  line  relaxation  schemes  in 
Chapter  IV,  the  high  frequency  (oscillatory)  components  of  the  error  are  extirpated 
much  more  efficiently  by  these  schemes  than  are  the  low  frequency  (smooth)  compo¬ 
nents.  The  result  of  applying  a  relaxation  scheme  to  generate  an  approximate  solution 
on  a  specific  grid  size  (call  it  a  fine  grid,  size  N),  is  that  the  oscillatory  components 
of  tlie  error  are  smoothed.  After  sufficient  work  has  been  done  on  the  fine  grid  to 
smooth  the  error  (in  effect,  when  the  relaxation  process  stalls),  the  problem  may  be 
shifted  to  a  coarser  grid  (grid  size  y)  where  the  smooth  components  of  the  error 
become  oscillatory  (see  figure  24,  taken  from  [Ref.  3]).  The  relaxation  scheme  is  then 
applied  again  to  smooth  the  oscillatory  components  of  the  error.  The  information 
gained  from  smoothing  on  the  coarse  grid  is  transferred  back  to  the  fine  grid,  where  it 
becomes  a  correction  to  the  original  approximation.  That  is,  the  result  of  smoothing 
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on  the  fine  grid,  transferring  to  the  coarse  grid,  smoothing  on  the  coarse  grid,  and 
transferring  back  to  the  fine  grid  is  to  produce  a  coarse  grid  correction.  This  pro- 
(•{'ss  of  using  multiple  grids  to  obtain  an  approximate  solution  has  the  advantages  of 
a.)  more  effectively  targeting  the  error  components,  and  b)  requiring  less  work,  since 
the  coarse  grid  has  only  2“*^  the  number  of  points  on  the  fine  grid,  where  d  is  the 
dimension  of  the  domain.  Multigrid  also  allows  for  ways  to  improve  upon  the  initial 
guess  that  is  used  in  the  smoothing  (see  [Ref.  3]).  However,  since  we  are  solving  a 
t  ime-stepping  problem,  where  the  current  solution  becomes  the  initial  guess  for  the 
next  time  step,  this  will  not  be  a  concern  for  us.  Moreover,  while  there  is  a  good  deal 
iiKue  to  multigrid  than  coarse  grid  correction,  this  concept  forms  the  foundation  of 
our  solution  process. 


k  =  4  wave  oq  N  =  16 


H 

a 

k  s!  4  wave  oq  N  =  8 

Figure  24.  A  wave  with  wave  number  =  4  on  {N  =  16)  is  projected  onto 
(iV  =  8).  For  N  =  16,  k  =  4  implies  that  the  wave  is  |  the  way  up  the  spectrum, 
while,  for  =  8,  A;  =  4  is  the  wave  that  is  |  =  |  way  up  the  spectrum.  Thus,  the 
coarse  grid  “sees”  a  wave  which  is  more  oscillatory. 

In  order  to  make  effective  use  of  this  technique,  we  must  specify  what  informa¬ 
tion  to  transfer,  and  how  to  transfer  it;  first  we  consider  what  information  to  transfer. 
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lve('all  that  we  are  attempting  to  solve  Equation  IV. 3, 

Mf-/, 

iiiid  siij^pose  that  we  have  an  approximate  solution  v\  the  error  is  given  by  e  =  T*  —  u , 
wliere  T*  is  the  exact  solution.  Thus,  the  error  satisfies 

Me  =  f  —  Mu 

-  r,  (VI.l) 

where  7^  is  the  residual.  The  smooth  components  of  the  error  are  the  troublesome 
ones,  since  the  relaxation  schemes  ’'kill”  the  oscillatory  modes.  Since  smooth  modes 
on  a  line  grid  become  oscillatory  when  projected  onto  a  coarse  grid,  it  is  natural  to 
( onsider  transterring  a  representation  of  the  error,  i.e.,  the  residual,  from  one  grid  to 
the  iiexth  In  this  way,  the  relaxation  schemes  used  on  the  coarse  grid  will  reduce 
components  of  the  error  that  could  not  be  reduced  on  the  fine  grid.  Additionally,  we 
know  that  relaxation  smooths  the  error,  and  therefore  we  can  accurately  represent 
the  error  on  the  coarse  grid.  After  transferring  the  residual  to  the  coarse  grid,  we 
can  relax  on  the  coarse  grid  version  of  the  residual  equation  (Equation  VI.  1),  solving 
for  the  error.  Thus,  when  we  have  determined  an  approximation  to  the  error  on  the 
coarse  grid,  we  can  transfer  it  back  to  the  fine  grid  as  a  correction  to  the  current 
approximation.  That  is,  since  T*  =  u  +  e,  if  we  know  v  we  can  correct  it  by  adding 
an  approximation  of  e.  This  process  can  be  outlined  as  in  the  following  steps,  where 
li  represents  the  grid  spacing  on  the  fine  grid,  H  the  spacing  on  the  coarse  grid,  and 
we  assume  2h  =  H . 

Coarse  grid  correction  (two-grid  scheme) 

•  Relax  on  M^T^  =  on  to  obtain  an  approximation  v^. 

^  There  are  a  number  of  other  good  reasons  to  transfer  the  residual,  as  outlined  in  [Ref.  3]. 
However,  there  are  other  multigrid  schemes  that  do  not  transfer  the  residual,  such  as  the  Full 
Approximation  Scheme  (FAS)  (see  [Ref.  2]  and  [Ref.  1]).  FAS  is  useful  in  dealing  with  nonlinear 
problems  but,  since  our  problem  is  linear,  we  will  not  employ  it. 
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•  Compute  the  residual  7^‘  = 

•  Solve  the  residual  equation  =  7“^  on  to  obtain  an  approximation 

to  the  error  . 

•  Correct  the  approximation  obtained  on  with  the  error  estimate  obtained 

on  +  e'^-lRef.  3] 

The  task  of  solving  on  can,  thus,  be  exchanged  for  the  task  of  relaxing  on 
and  then  solving  on  0,^ .  The  requirement  to  solve  on  can  be  treated  in  like 
fashion:  the  task  of  solving  can  be  exchanged  for  the  task  of  relaxing,  and  then  solving 
(ui  (.he  next  coarser  grid.  This  procdeure  of  transferring  to  successively  coarser  grids 
can  be  continued  until  a  grid  is  reached  011  which  an  “exact”  solution  is  possible. 
That  is.  the  solution  can  be  generated  by  recursive  use  of  the  coarse  grid  correction, 
which  can  be  represented  as  follows: 

Solve  on  by  relaxing  oii  and  solving  on 

Solve  on  by  relaxing  on  and  solving  on 

Solve  on  by  relaxing  on  and  solving  on  1)*^. 

“Exact”  solve  on  coarsest  grid. 

Correct  on  0“*^. 

Correct  on 
Correct  on 

This  process  of  starting  on  a  fine  grid,  moving  to  the  coarsest  grid,  and  then  going 
Imck  to  the  fine  grid  is  known  as  a  V-cycle  (see  Figure  25).  Each  time  the  transfer 
to  a  coarser  grid  is  made,  the  relaxation  process  there  smooths  another  portion  of 
the  error  component  spectrum  (see  Figure  26).  Thus,  by  the  time  the  coarsest  grid 
is  reached  and  the  problem  has  been  solved,  all  of  the  error  components  have  been 
smoothed.  When  the  problem  is  solved  on  the  coarsest  grid,  the  approximate  error 
is  transferred  to  the  next  finer  grid  to  become  a  correction.  The  corrected  error  is 
then  passed  to  the  next  finer  grid,  and  so  on  until  we  are  back  on  the  finest  grid, 
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wliere  the  original  approximation  is  updated  with  the  composite  correction  from  the 
roarser  grids.  While  there  are  other  multigrid  methods  which  may  prove  more  useful. 
ile|)('nding  on  the  nature  of  the  problem  (see  [Ref.  3]),  we  use  the  V-cycle. 


Figure  25.  Schedule  of  grids  for  a  V-cycle. 
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Figure  26.  Error  component  spectrum  for  various  grid  sizes.  On  the  fine  grid,  /i, 
relaxation  smooths  the  high  frequency  components  of  the  error  (heavy  black  line). 
When  the  problem  is  transferred  to  the  next  coarser  grid,  2/i,  the  high  frequencies 
of  the  remaining  (unsmoothed)  components  (heavy  black  line)  are  smoothed.  This 
process  continues  all  the  way  down  to  the  coarsest  grid,  where  the  problem  is  solved, 
eliminating  the  remaining  components  of  the  error.  By  continuing  to  transfer  to 
successively  coarser  grids,  all  frequency  components  of  the  error  are  smoothed. 

So  far,  we  have  not  indicated  how  information  will  transferred,  nor  how  an 
error  estimate  on  the  coarse  grid  (solving  will  be  computed.  We  will 
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(leal  with  the  former  question  in  the  next  section;  the  latter  question  revolves  arouini 
(lecitling  what  to  use  as  the  coarse  grid  operator.  M^.  One  approach  is  to  simply 
discretize  the  problem  on  the  coarse  grid  in  precisely  the  same  fashion  that  it  is 
discretized  on  the  fine  grid.  This  has  the  advantages  that  the  work  has  already  been 
done,  and  that  it  is  easy  to  implement  computationally.  While  this  method  generally 
woi  k.s  (on  model  problems),  it  has  the  disadvantage  that  it  is  often  not  true  to  the 
physics,  e.g.,  conservation  may  no  longer  be  enforced.  There  are  other  methods  for 
determining  the  coarse  grid  operator  which  are  true  to  the  physics  and  allow  for  strong 
theoretical  treatment  of  convergence  and  other  properties  (see  [Ref.  2]  and  [Ref.  3]). 
However,  because  discretization  on  the  coarse  grid  generally  works,  and  because  it 
simplifies  the  computations,  forming  in  this  manner  is  the  method  that  we  use. 

We  conclude  this  section  with  a  discussion  of  what  happens  to  boundary  con¬ 
ditions  when  we  transfer  the  residual.  The  boundary  conditions  for  our  problem  are 
a  mixture  of  Dirichlet  and  Neumann  conditions  (Equations  11.18,  11.19,  and  11.20). 
By  relaxing  on  Equation  IV. 2, 

L{T]  =  /, 

(where  we  have  dropped  the  superscript  h)  we  obtain  an  approximation  u,  which  gives 
rise  to  the  residual  equation,  r  —  f  —  L[u].  Suppose  that  T*  is  the  exact  solution,  so 
that  /  =  L[T*]]  we  require  that  the  approximation  satisfy  the  same  conditions  as  the 
exact  solution  on  the  boundary,  d^l.  Thus,  the  residual  equation  becomes 


r  =  L[r]-L{v] 

=  [ -  aKV^)rdV]  -  [ /^(i  - 

=  i  -  v)dV  -  UK  -  V'‘v)dV 

=  -  j^(r -v)dV -aK  j^(V'r -Vv)-ndS. 


(VI.2) 


Therefore,  since  the  boundary  conditions  satisfy  (T*  =  u)|an  and  all 


-  ^  =  V  ■  n  denotes  the  outward  normal  on  dii. 
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boundary  conditions  for  ?•  become  homogeneous;  the  first  term  on  the  right  hand  side 
of  Equation  VI. 2  indicates  that  this  is  true  for  the  Dirichlet  conditions,  the  second 
term  indicates  that  it  is  true  for  the  Neumann  conditions.  Thus,  since  we  are  solving 
rlie  lesidual  equation  on  all  the  successively  coarser  grids,  we  impose  homogeneous 
lioundary  conditions  on  all  but  the  finest  grid. 

A.  INTERGRID  TRANSFERS 

In  order  to  transfer  information  between  grids,  we  develop  operators  that 
restrict  data  from  a  fine  grid  to  a  coarse  grid,  and  interpolate  data  from  a  coarse 
grid  to  a  fine  grid  (see  [Ref.  3]).  These  operators  are  designated  by  IsubscHpt*^'’  where 
1  lie  suliscript  indicates  the  grid  from  which  the  information  is  being  transferred,  and 
t  he  superscript  indicates  the  grid  to  which  the  information  is  being  transferred.  Thus, 
the  restriction  operator  is  indicated  by  ij/,  since  the  information  goes  from  the  fine 
grid  to  the  coarse  grid,  and  the  interpolation  operator  is  represented  by  I^,  since 
the  information  goes  from  the  coarse  grid  to  the  fine  grid.  In  other  words,  restriction 
is  a  process  of  determining  what  values  to  assign  to  coarse  grid  points,  based  on  the 
values  at  fine  grid  points;  interpolation  is  a  process  of  determining  what  values  to 
assign  to  fine  grid  points,  based  on  the  values  at  coarse  grid  points.  In  [Ref.  1],  such 
o|)erators  are  developed  which  take  advantage  of  FVE  characteristics.  We  consider 
two  types  of  restriction  operators,  one  which  follows  from  the  physics  and  one  which 
is  very  simple,  and  one  type  of  interpolation  operator. 

One  of  the  advantages  of  the  FVE  method  is  its  fidelity  to  conservation  consid¬ 
erations.  The  restriction  technique  that  we  describe  first  follows  from  the  conservation 
notion.  In  order  to  determine  what  value  to  assign  to  a  coarse  grid  point,  the  fine  grid 
is  laid  over  the  coarse  grid,  where  the  coarse  grid  control  volumes  align  with  fine  grid 
lines  (see  Figure  27).  The  control  volume  for  a  given  coarse  grid  point  includes  all  or 
part  of  the  control  volumes  of  nine  fine  grid  points.  Hence,  the  value  of  a  quantity  on 
the  coarse  grid  includes  contributions  from  the  value  of  the  quantity  on  each  of  the 
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nine  line  grid  points.  For  each  of  these  fine  grid  points,  the  contribution  is  determined 
by  computing  the  fraction  of  the  associated  fine  grid  volume  that  is  contained  in  the 
coarse  grid  volume.  The  value  of  the  quantity  at  the  fine  grid  point  is  multiplied  by 
this  fraction,  and  the  result  is  the  contribution  to  the  coarse-grid  value.  The  value 
i)f  the  quantity  at  the  coarse  grid  point  is  thus  the  sum  of  the  nine  contributions 
determined  in  this  manner. 


In  the  case  of  cylindrical  coordinates,  care  must  be  taken  when  deciding  what 
percentage  of  the  fine  grid  volumes  fall  within  each  coarse  grid  volume.  Since  the 
volumes  change  with  the  radius,  the  restriction  operator  will  be  a  function  of  grid 
location.  The  stencil  for  the  restriction  operator  gives  the  percentage  of  the  fine  grid 
value  that  is  assigned  to  the  coarse  grid  value  based  on  the  percentage  of  the  fine 
grid  volume  that  falls  within  the  coarse  grid  volume  (see  Figure  28).  For  example, 
the  amount  of  the  value  of  the  point  at  (k,j  +  1)  that  is  used  in  the  sum  is  found  by 
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Figure  28.  The  fine  grid  region,  comprising  contributions  from  the  nine  fine  grid 
control  volumes,  that  corresponds  to  a  coarse  grid  control  volume.  The  central  fine 
grid  point  is  labeled  {k,j  even)  for  which  the  radial  distance  r  =  j7i,  and 

corresponds  to  coarse  grid  point 


first  determining  how  much  of  its  fine  grid  control  volume  (right-center,  or  rc)  falls 
within  the  region  that  corresponds  to  the  coarse  grid  control  volume.  This  volume 
is  then  divided  by  the  total  control  volume  for  the  point  {k,j  -F  1)  to  determine  the 
volume  percentage.  So,  if  the  fine  grid  location  is  at  the  radial  distance  r  =  jh,  then 
we  have  for  the  volume  7'c  (a  toroidal  prism) 

^hl{j  +  i)V-U  +  \Th% 

or 

which  is  divided  by  the  fine  grid  control  volume  for  the  point  at  {k,j  -f  1), 

+  -  U  + 

or 

Trh^{2j  -F  2), 
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giving 

J  +  i 

2{J  +  1)' 

Similarly,  the  ratio  for  the  left-center  (Ic)  portion  of  the  stencil  is 

J-1 

20’ -D' 

and  the  remainder  of  the  stencil  is  given  by, 

'i(^)  I  i{^)' 

irld)  i  i(Z±i) 

_  4V  j-1  /  2  4  V  j-i-1  J 

This  is  a  stencil  for  interior  points;  similar  calculations  determine  the  stencils  for 
hoiinclary  points. 
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Figure  29.  Neighboring  coarse  grid  points. 


Since  this  operator  is  both  somewhat  out  of  the  ordinary  and  dependent  on  grid 
location,  we  consider  it  worthwhile  to  try  to  get  a  sense  of  how  this  restriction  operator 
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Heals  the  quantities  being  restricted.  For  this  purpose,  we  consider  a  comparison 
between  the  values  assigned  to  neighboring  coarse  grid  points,  as  in  Figure  29,  when 
the  restriction  operator  acts  on  a  grid  with  constant  value  Tk.j  =  1.  The  value  for 
I oarse  grid  point  number  1  is  given  by 


2  4-:^ 


+  2  + 


J  +  l 


7^  -  -  7^-14-^ 

2  +  2f - - - i - )  =  O  +  Pf-Ji _ LZjL— 


=  4  + 


The  corresponding  value  for  coarse  grid  point  number  2  is 

3  ■  “  2-1 


p-^j  +  T 


=  4  + - i - . 

2(i  -  3)(i  -  1) 

■Since  (j  +  1)  >  (j  —  3)  for  2  <  j  <  N  —  2,  then  which  impli 


and  therefore 


2(7  +  l)(i-l)  2(j-3)U-l)’ 


1  1 

‘‘+20  +  1)0-1)  ^'‘+20-3)0-1)' 


Thus,  if  the  restriction  operator  acts  on  a  constant  vector,  the  values  assigned  to  the 
coarse  grid  points  decrease  as  the  radial  distance  increases. 

Does  this  make  sense  in  terms  of  the  physics  of  the  problem?  Should  a  constant 
function  remain  a  constant  after  restriction,  or  not?  In  general,  one  might  expect 
that  restriction  would  transfer  a  constant  to  a  constant.  The  conservaton  restriction 
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o|)(Tat()r  acts  on  the  basis  of  volume  percentages  and,  thus,  considers  the  values 
assigned  to  each  control  volume  in  terms  of  the  actual  volume.  In  some  sense,  this 
restriction  operator  converts  the  values  assigned  to  a  control  volume  into  a  “density”, 
and  tlien  transfers  this  amount.  In  the  case  of  cylindrical  coordinates,  assigning  a 
constant  value  to  each  control  volume  means  that  the  ratio  of  the  value  to  its  volume 
decreases  with  radial  distance.  Thus,  the  result  of  the  above  computation  is  no 
surprise.  Moreover,  it  seems  to  make  sense  physically  to  transfer  “densities”  in  terms 
of  enforcing  conservation.  However,  in  order  to  regain  the  fine-grid  interpretation  of 
the  information  after  it  has  been  transferred  may  require  more  work.  In  particular, 
some  compensation  may  need  to  be  made  relative  to  the  coarse  grid  volumes  (i.e., 
convert  the  “density”  back  to  the  original  type  of  information).  The  question  of  how 
of  whether  to  make  this  compensation  will  not  be  addressed  here.  Additionally,  it  is 
certainly  possible  that  the  information  to  be  transferred  does  not  make  sense  in  the 
context  of  densities.  In  that  case,  physical  intuition  may  need  to  be  suspended  while 
the  efficacy  of  the  solution  process  is  determined. 

The  other  restriction  operator  we  introduce  is  the  injection  operator.  Injection 
is  accomplished  by  simply  a.ssigning  values  to  the  coarse  grid  points  directly  from  the 
corresponding  fine  grid  points,  vffj  =  v^kaj  [Ref.  3]).  The  injection  operator 
has  the  advantage  of  being  a  linear  operator,  whereas  the  conservation  restriction 
operator  may  not  be  linear.  Note  that,  when  the  conservation  restriction  operator  is 
applied  to  a  grid  with  constant  value  Tkj  =  1,  the  result  is  T2k,2j  ~  4.  If  the  injection 
restriction  operator  is  applied  to  the  same  grid,  T2k.2j  =  1-  Thus,  care  must  be  taken 
when  comparing  these  two  restriction  operators. 

In  order  to  transfer  from  the  coarse  grid  to  the  fine  grid,  we  must  choose  an 
interpolation  operator.  Since  we  are  assuming  that  the  basis  functions  for  the  FVE 
method  are  linear,  a  natural  choice  is  linear  interpolation.  A  continuous  -piecewise 
linear  function  is  also  a  continuous  /j-piecewise  linear  function.  That  is,  the  “coarse” 
finite  element  space  is  a  subspace  of  the  “fine”  finite  element  space  (see  [Ref.  1]). 
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rilis  means  that  linear  interpolation  corresponds  to  the  “natural  embedding"  from 
the  coarse  finite  element  space  to  the  fine  element  space.  Its  advantages  are  its 
simplicity  and  its  linearity.  The  stencil  is  given  by 


B.  AMPLIFICATION  MATRIX 

.Now  that  we  have  our  fine  grid  operators,  coarse  grid  operators,  and  intergrid 
transfer  operators,  we  conduct  analysis  using  the  DFT  (see  Appendix  A)  similar  to 
that  in  Chapter  IV  to  obtain  some  indication  as  to  how  a  two-level  multigrid  scheme 
using  these  operators  will  perform.  We  start  by  noting  that  the  process  of  relaxing 
1)11  the  fine  grid,  obtaining  a  coarse  grid  correction,  and  relaxing  again  on  the  fine 
grid  can  be  represented  by  the  two-level  error  propagation  matrix,  P'"'^ {C G) P''' .  Note 
that,  while  h  as  a  superscript  identifies  grid  spacing,  i/j  and  1^2  are  exponents.  Thus, 
P"'  and  P''^  represent  relaxation  and  U2  times  on  the  fine  grid,  and  CG  Is  the 
coarse  grid  correction.  Recall  that  coarse  grid  correction  is  the  following  series  of 
steps: 

Coarse  grid  correction  (two-grid  scheme) 

•  Relax  on  on  to  obtain  an  approximation 

•  Compute  the  residual  7’^  =  /^  — 

•  Restrict  the  residual 

•  Solve  the  residual  equation  on  0^  to  obtain  an  approximation  to 

the  error  . 

•  Interpolate  the  error  . 

•  Correct  the  approximation  obtained  on  with  the  error  estimate  obtained 

on  0^:  4—  -b  e^. 
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rims,  CG  can  be  written  as 


CG  =  I  - 

where  /  is  the  identity  operator,  ij^  and  are  the  restriction  and  interpolation 
operators,  is  the  fine  grid  operator,  and  is  the  coarse  grid  smoother.  Each 

of  the  operators  in  the  two-level  error  propagation  matrix  is  expanded  in  a  DFT 
(see  .Appendix  A  and  Chapter  IV),  making  the  resulting  expansions  functions  of 
l.in  =  — Y  +  I  :  y,  the  indices  of  the  expansion  in  the  frequency  domain.  These 
expansions  are  used  to  construct  an  amplification  matrix,  A{l,m.),  which  is  also  a 
function  of  I  and  m.  The  largest  spectral  radius  of  this  matrix  over  the  values  that 
l.  iu.  have  on  the  coarse  grid,  i.e.,  /,m  =  —  ^  -h  1  :  y  (this  relationship  is  discussed 
later),  will  give  us  a  two-level  asymptotic  convergence  factor,  A.  This  factor, 
much  like  the  values  determined  for  the  relaxation  schemes,  will  give  a  bound  on 
how  well  the  two-level  scheme  can  be  expected  to  perform.  In  order  to  guarantee 
convergence,  A  <  1  is  necessary;  the  smaller  A  is,  the  better. 

In  Chapter  IV,  we  used  DFT  expansions  of  the  error  equations  derived  from  the 
relaxation  operators  to  determine  the  frequency  domain  coefficients  for  the  operators. 
VVe  now  use  the  same  technique,  again  using  error  equations,  to  determine  the  matrix 
entries  for  the  operators  in  the  two-level  scheme.  The  amplification  matrix  is  formed 
by  multiplying  together  these  operator  matrices, 

A(/,  m)  =  (P'‘)‘^^(I  -  ,  (VI.3) 


where 

•  is  the  4x4  matrix  for  the  fine  grid  discretization  operator. 

•  P^  is  the  4x4  matrix  for  the  fine  grid  relaxation  operator. 

•  (L^)”'  is  the  1x1  matrix  for  the  coarse  grid  relaxation  operator. 

•  is  the  1x4  matrix  for  the  restriction  operator. 

•  is  the  4x1  matrix  for  the  interpolation  operator. 
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•  I  is  the  4x4  identity  matrix,  (see  [Ref.  2]) 

We  go  through  the  details  of  computing  the  matrix  entries  for  the  fine  grid  operator, 
results  for  the  remaining  operators  are  then  presented. 

From  the  FVE  stencils  for  U’’,  we  have  the  operator  equation 

k(new)  _  r  h  k(old) 

^  2k:2j  —  ^  ^’2k:2j  ’ 


^  (  (64i  -  5)t>2;“l2,  ' 

^  +448^41:'''  +(«4j  +  11)41*1.. 

+(.12j  —  ■l-(64^  +  ■'1)41'1i!2j  j 

(  4>ri22  'l 

+(-i+4,>ssl.  -16.41:1?  +(i+4.')41:iii,  . 

\  '^J^2k-l,2j  / 

where  the  superscript  li  indicates  a  vector  component  on  the  fine  grid,  {old)  and 
(new)  indicate  before  and  after  the  operator  acts,  and  subscripts  2k  and  2j  are  the 
respective  fine  grid  row  and  column  indices.  These  indicial  values  are  used  because 
we  shortly  will  need  to  discuss  the  correspondence  between  the  fine  and  coarse  grid 
points;  only  those  fine  grid  points  with  indices  2k  and  2j  correspond  to  coarse  grid 
points  (with  indices  j  and  k).  Ordinarily  (see  [Ref.  2]),  the  consideration  of  indicial 
values  is  not  a  concern.  In  this  case,  however,  since  the  stencils  are  grid  location 
dependent,  and  since  the  analysis  of  the  amplification  matrix  must  apply  to  both 
the  fine  and  coarse  grids,  the  j  that  is  used  is  the  column  index  on  the  coarse  grid, 
requiring  that  2j  be  the  column  index  on  the  fine  grid.  Thus,  expanding  both  sides 
in  a  DFT,  we  have 
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When  we  transfer  information  from  the  fine  grid  to  the  coarse  grid,  we  must 
remember  that  not  all  of  the  freqencies  on  the  fine  grid  can  be  represented  on  the 
('oai  se  grid.  The  highest  frequency  that  can  be  resolved  on  any  grid  is  the  Nyquist 
frequency  (see  [Ref.  12]),  /  =  2^,  where  Ax  is  the  grid  spacing.  Hence,  on  the 
coarse  grid,  ,  Ax  =  ^  and  f  =  Therefore,  in  order  to  make  this  analysis 
germane  to  both  fine  and  coarse  grids,  we  must  rewrite  the  sums  in  the  expansions 
so  that  |/|,  |mj  <  This  can  be  done  as  follows: 


N  N 


K  K 

4  4 


X:  E  H,„c(i.m)=  y;  (HmCc/.m)  (vi.4) 

N  N  N  N  ^ 

m  +  — )  +  +  y)  • 


76 


Thus,  rewriting  the  DFT  e,\pansion  of  the  operator  equation  to  consider  those 
rr('(|iu'nc.ies  that  can  be  represented  on  the  coarse  grid  and,  making  use  of  the  or- 
1  liou.unality  i)roperty  of  the  complex  exponential  (see  Appendix  A),  we  can  equate 
individual  terms  of  the  sum.  a.nd  then  divide  by  C{l,m)  7^  0  to  give 

vS!  =  [Jj((64,-5)V£.^+(32,+5)V'£e^ 

+(64i  -  n)u':;c=^  +448iv;!:; + (64j  + 

+(.'«i  -  5)v;1.T-^  +  (644  + 

-^(44V-i‘::!,‘^  +  (-l+44)VS,-‘^ 


-164VS!+(l+44)V£.lc‘4S«+4,Vjy=f= 


or,  factoring  out  V; 


/  i2nl  i27r(t  +  m) 

=  ^  -  5)e—  +  (32i  +  5)e^^ 

+(64j  —  ll)e  N  +  448j  +  (64j  +  ll)e" 

t27r(f  +  m)  ^  i2irl  \ 

+(32j;  —  5)e  ^  +  (64j  +  5)e  ^  1 

aKh ,  .  t2ffi  i27rm 

- ^(4je  w  4.(_i+4j)e  -v 

-16i  +  (1  +  4i)e'^  +  4je“^)]  V,^^, 


yin) 


f  ,  .  .2>r(  ,  .  ,  ■2ir((4n 

[-(64i  -  5)e  N  -  (32j  +  5)e  n 
+(64j  —  ll)e  w”*  +  448j  +  (64j  +  ll)e“ 


-(32i  -  5)e- 


»2ir((+m) 

N 


(64j  +  5)e 


OiKh,  (  ,  .  i2nl  ,  , 

- —  (-4je  w  +  (-1  +  4j)e  ^ 

-16j  +  (1  +  4j)e^  -  4je“^)]  ,, 
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since 

ilEL  iT,  i2JlL  ,  i2^(m4-4)  ,2trm 

e  ^  =  e  ^  e  =  —e  ^  ,  and  e  ^  =  — e  ^  . 

The  matrix  of  coefficients  for  the  fine  grid  operator  is  a  4x4  diagonal  matrix, 
(see  [Ref.  2]).  To  see  why  the  matrix  is  4x4,  recall  that  the  sum  in  the 
original  DFT  expansion  of  the  operator  equation  was  rewritten  so  that  the  indices 
satisfy  |/|,  |m|  <  One  result  of  this  is  to  rewrite  each  of  the  individual  terms  of 
the  sum  as  a  combination  four  terms  (Equation  VI. 4).  That  is,  for  each  dimension 
of  the  problem,  represented  by  I  and  m,  there  are  two  components  of  each  individual 
term  in  the  sum,  which  are  indexed  by  |/|,  |m|  <  ^  and  ^  <  |/|,  \m\  <  y.  Thus,  the 
matrix  of  coefficients  is  a  2^x2^  matrix,  where  each  row  in  the  matrix,  say  the  ith  row, 
comprises  the  coefficients  used  in  determining  the  corresponding  ith  component  in  the 
new  vector  as  a  linear  combination  of  components  of  the  old  vector.  In  this  case  the 
matrix  is  diagonal,  since  for  each  of  the  four  components  of  the  new  vector,  the  only 
contribution  from  the  old  vector  comes  from  the  corresponding  component.  In  other 
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,  i2ni  »2)r(l+m) 

A3.3  =  ^  (^(64;  -  5)e—  -  (32i  +  5)e-^^ 

-(64j  -  ll)eT  +  448i  -  (32;  +  ll)e^ 

t23r(i+m)  i7wi  \ 

—(32;  —  5)e“  ^  +  (64;  +  5)e  ^  j 

anh  /.  .  /  -  .  .X  »2irm 

— —[4jeN  _(_i+4j)e-  n 

-16;-(l+2;)e  ^  +  4;e  n  jj  ^ 
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A4,4  =  [-(64j  -  o)e  N  +  (32j  +  D)e  ^ 

-(64i-  ll)e^+448i-(64i  +  ll)e'^ 

/  ^  .  \  t27r(/-f  m)  ^  x2'kL  \ 

+(32j  —  o)e  'V  _  4-  5)e  n  j 

a«;/i  /  .  *27r<  t27rTn 

- j-  (^-4je  -  (-1  +4j/)e  n 

.  i27rm  — i27r^  \  ' 

-16j-(l+44)e  -v  -4je  j  . 

A  similar  process  must  be  applied  to  determine  each  of  the  other  matrices 
that  comprise  A{l,m.).  In  order  to  compute  the  entries  for  the  matrix  P^\  recall  from 
(diapter  IV,  Equation  IV. 25,  that  vve  can  write  the  error  equation,  =  Pe°^‘\  as 


OckIx  pieu;)  _ 


f  (new)  .  (Tiett;)  7-)  (new)  ^ 

^^2k-2j^l  ‘2j  ^^2k-\-\:2j+l) 

2  V  ^^2fc-l,2i  —  ”'2Jk,2j-l 
^^2A:,2j+l  ^^2fc+l,2i;) 

where  .4  =  32;  -5,5  =  64;  +  5,  C  =  64i  -  11,5  =  448;,  5  =  64i  +  11,5  = 
fi4ji  —  5,  G  =  32ji  +  5,  and  5  =  4j,  /  =  —  1  +  4j,  J  =  —  16j,  /<'  =  !+  4j,  L  =  44.  Using 
similar  calculations  to  those  above,  we  have  the  four  entries,  IT, of  the  diagonal 
matrix  P^‘: 


n,.i  = 


n.>  2  = 


r  T-,  ^27rm  t2nl  „  -■  ,vk;/i  r  r-  r  1^1 

1  — £^e  N  ^  Fe  N  -  Ge  n  J  -  Ae  ^  -Lcn 

^  i2nl  ^  —i2irm  _ u  r  — t27r£  ,  »27rm  ,,  ’ 

+  5e"^  +  Ce—^  +  D]-  ^[He-sr-  +  /e“^  +  J] 


i2‘K{l  +  m 


ith! 


r  .  »27r(Hm) 

:[Ae  N 


r  1-1  *27rm  —  i2nl  ^  « 

[— 5e  w  +  Fe  N  +  (?e' 


tVK/ir  r/-  »27rTn  ^  i2irl . 

—  ^[— A  e  w  +  Le  N 


^  i2Tr(i+»Ti)  ^  i27ri  ^  --«27rm  _ i.  ,  ,,  — i^irZ  ,  i27rm  ’ 

-Ae^  ^  —  Be  +  Ce  N  + n  +  le  ^  +  J] 


n3,3  -- 


y  3  r  ^ 
i2ir(Z+m) 


77,  ilni  ^  t27r(Hm)  aKhir^  »24rm  ^  i27ri^ 

Fe  N  +  Ge  n  j  ^  -^[Ae  n  ^Lcn 


j-3  P  i2ir(Z+m)  i27T-i  ^  -i2irm  ^  -*27ri  ^  t27rm 

^[-Ae - N  +  Be  —  -  Ce^r-  +  5]  -  ^[5e-^  -  /e  ^  +  J\ 


114.4  — 


l3  <27rm  _  t27rf  ^  t27r(l+m)  ,  »2Trm  - 

^Fe  N  Fe  N  -  Ge  n  ]  -  2^[Fe^r-  +  Le 


_  i2ir<  ^  — i27rm  _ tr  --  — •2'7rl  ,  i2jrm  ...-i 

Be-—  -  Ce-ir-  +  D]  -  ^[-He—  -  /e"—  +  J] 
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In  similar  fashion,  the  entry  for  the  1x1  matrix  (L^)  \  corresponding  to  the 
coarse  grid  relaxation  operator,  is  given  by 


rrU  »27rm  ^  tt  t2nt 

-  -^[-A^e—  -  L^e— 


^lA^e- 

■f.itiA  e 


i27t(  i  +  m) 

N 


+  BHe-^  +  +  D"]  -  +  ./"]' 


.’here  =  iej-5,  5"  =  32i  +  .5.  =  32i-  11,  D"  =  224;,  =  32;  +  11,  = 


32;  -  5.  =  16;  +  5,  and  =  2;,  /^  =  - 1  +  2;,  =  -8;,  =  1  +  2;,  =  2;. 

The  equation  for  the  conservation  restriction  operator,  /^,  is 


= 


_  4)  f,./t  1  _L  0+  I 

“  4^2;  —  1)  ^^2k,2i-\  ■r  ^2/c+l,2i-l/ 

{^2fc-l,2i  "b  ‘^'^2k,2j  "b  ^2*;+l,2i} 


(2j  +  f 


2fc,2i+l  ^  ^2A:+l,2j+l 


Let  =  ■(2JZ1)  ^  “  W+if  expand  this  relation  in  a  DFT.  Considering 

those  frequencies  that  can  be  represented  on  the  coarse  grid  and,  making  use  of  the 
orthogonality  property,  we  can  equate  individual  terms  of  the  sum.  We  divide  by 
C  (<,m)  =  e  N  e  N  ,  the  coarse  grid  analog  to  C{l,m),  to  give  the  components  of 
the  matrix  of  coefficients  corresponding  to  the  restriction  operator,  =  (Ai,,).  This 
matrix  is  1x4  since  the  output  of  the  operator,  a  single  component  on  the  coarse  grid, 
is  a  linear  combination  of  2^  components  on  the  fine  grid.  The  entries  in  the  matrix 
are  as  follows: 

A  ^  ~b  cos  (  pj  i  r  i2irm  i2Km  1 

Ai,i  =  - [ge-— +  2  +  Ae— j  , 

A  ^  i2trm  j2»rm  1 

Ai,2  =  - [ge-— +  2  +  Fe— J  , 

A  ^  )  r  y-k  *27rm  ^  _  t27rm  1  , 

Ai,3  =  - +2- Re-^\,  and 

-  1-  cos  { ip  i2nm  i2mn  I 

Ai,4  =  - Y^[-Qe-— A-2-Re—\. 
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The  equation  for  the  injection  restriction  operator  is  —  vl^f^  2j  which,  using 
similar  calculation,  results  in 

V*  =  V'l  +  ^ 

so  that  =  [1  1  1  1]. 

The  linear  interpolation  operator  gives  rise  to  the  following  relations: 


7/‘ 


^2A:+l,2i  “ 


■^2^^2J  +  l 


^2A:+l,2jf  +  l 


^k,j  + 


+1 


H  H 
'^k,j  I 


Expanding  these  in  a  DFT,  considering  the  coarse  grid  frequencies,  and  making  use 
of  the  orthogonality  property,  we  can  equate  terms  of  the  sums  and  divide  by  C'(/,  m) 
to  give 


rfi  _ +  e  ^  I 


— 2 — ^  =  v;*„  - 

/  »27rm  —i2nm  . 

\/H  ^  ^  ^  j  T//i  I  w/i  I  T/h  I  T//i 

/  >27rf-hm  ~t27r(f+m) 

/H  'S _ ^ _ !! _ i  _  t/a  I  T/^ ..  _L  lA  _ L  \A .. 


2  2 


The  components,  G^j,  of  the  matrix  of  coefficients,  I^,  are  found  by  solving  the  above 
system  of  four  equations  for  the  fine  grid  components,  such  that  =  The 

matrix  is  4x1  since  the  output  of  the  operator,  the  2^  components  of  the  fine  grid 
vector,  is  the  result  of  interpolation  operator  acting  on  the  single  component  of  the 
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(■i)ars('  grid  vector.  That  is, 


’  e,,,  ■ 

'  1  +  co«(^)  +  003(252)  +  cos(M!^)  ' 

02.1 

1 

1  -  003(^)003(252)  _  003(22152^) 

03,1 

~  4 

1  -L  cos(^)  -  cos{-^)  -L  cos(;^^) 

04,1 

_  1  -  003(f)  -  003(2221)  +  COSC-^)  ^ 

Thus,  we  have  all  the  elements  to  construct  the  amplification  matrix,  A{l,r7i), 
( hcjuation  VI. 3).  VVe  now  have  the  tools  to  analyze  the  expected  perfomance  of  a 
two-level  scheme. 

C.  NUMERICAL  RESULTS 

We  now  conduct  a  number  of  numerical  experiments  to  investigate  the  perfor¬ 
mance  of  the  multigrid  method,  again  using  Matlab.  First  we  examine  the  behavior 
of  the  amplification  matrix,  A(/,m),  and  the  largest  spectral  radius  of  the  family  of 
matrices  A(/,  m)  for  all  I,  m.  We  then  present  the  results  of  the  repeated  use  of  the 
V-cycle  to  solve  Equation  IV. 2,  as  compared  with  the  iterative  solutions. 

In  order  to  make  use  of  the  amplification  matrix,  we  must  choose  a  relaxation 
method,  and  the  number  of  times  to  relax  on  the  way  down  and  on  the  way  back  up, 
i.e.,  we  must  choose  i/j  and  1^2  in  Equation  VI. 3.  We  experiment  first  with  Gauss- 
Seidel  relaxation  and  three  different  combinations  of  (1^1,  ^'2)  V-cycles:  a  (0, 1)  cycle, 
a  (2,1)  cycle,  and  a  (10,10)  cycle  (recall  that  a  =  1,  meaning  fully  implicit  time 
stepping,  and  the  time  step  g  =  h).  The  first  two  types  of  cycles  are  fairly  common 
(see  [Ref.  3]),  while  the  latter  is  included  for  reference.  The  information  presented  in 
Tables  VIII,  IX,  and  X  indicates  the  asymptotic  convergence  factor,  A,  for  specific  grid 
sizes,  grid  positions,  and  type  of  restriction  operator,  either  conservation  or  injection. 
Linear  interpolation  is  used  as  the  interpolation  operator  in  all  of  these  experiments. 
Tables  XI  and  XII  indicate  the  asymptotic  convergence  factors  for  horizontal  and 
vertical  line  relaxation  and  a  (2, 1)  cycle. 

Table  VIII  indicates  the  asymptotic  convergence  factors  that  result  from  a 
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(0,1)  cycle,  Gauss-Seidel 

Conservation 

Injection 

N  =  16 

N  =  32 

iV  =  64 

N  =  16 

N  =  32 

N  =  64 

2 

.9174 

.9651 

.9843 

.9397 

.9710 

.9858 

.8308 

.9111 

.9544 

.9144 

.9553 

.9771 

N  -  2 

.7472 

.8595 

.9259 

.8929 

.9422 

.9700 

■fable  VIII.  A  symptotic  Convergence  Factors  for  the  (0,  1)  Cycle  Using  Gauss-Seidel 
Relaxation. 

single  iteration  sweep  on  the  fine  grid  coupled  with  the  course  grid  correction  for 
l)oth  the  conservation  and  injection  restriction  operators.  In  general,  the  conservation 
r(‘striction  operator  seems  to  indicate  better  performance,  since  its  convergence  factors 
ai  e  lower  than  those  for  the  injection  operator  (recall  that  the  asymptotic  convergence 
factor  is  the  largest  spectral  radius  of  the  family  of  amplification  matrices,  which 
comes  from  the  DFT  expansion  of  each  component  in  the  two-level  error  propagation 
matrix).  Similar  to  the  analysis  on  relaxation  schemes  in  Chapter  IV,  the  convergence 
factors  decrea.se  with  radial  distance,  indicating  an  improvement  in  performance,  but 
increase  with  grid  size.  Even  though  the  convergence  factors  for  all  three  grid  sizes  are 
less  than  one,  indicating  that  convergence  is  very  likely  (but  not  necessarily  certain, 
since  these  factors  tell  us  nothing  about  what  is  going  on  at  the  boundaries),  at  least 
some  of  the  factors  for  all  three  grids  are  relatively  high,  indicating  that  convergence 
may  be  quite  slow  with  the  (0, 1)  cycle.  A  point  of  interest  is  to  compare  these 
results  with  those  in  Table  IV,  which  indicate  the  maximum  ratio  of  new  to  old  error 
components  for  Gauss-Seidel  relaxation.  One  might  expect  the  results  in  Table  VIII 
to  be  universally  better  than  those  in  Table  IV,  however  this  is  not  the  case.  In 
fact,  on  a  grid  of  size  N  =  32,  the  largest  asymptotic  convergence  factors  for  both 
the  conservation  and  the  injection  restriction  operators  are  larger  than  the  largest 
ratio  for  just  a  single  Gauss-Seidel  sweep  on  the  fine  grid.  In  other  words,  a  single 
Gauss-Seidel  relaxation  sweep  is  predicted  to  do  better  than  a  (0, 1)  cycle.  This  is  a 
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(2,1)  cycle,  Gauss-Seidel 

Conservation 

Injection 

2i 

N  =  16 

N  =  32 

iV  =  64 

yv  =  16 

N  =  32 

Af  =  64 

2 

.8223 

.9133 

.9574 

.8423 

.9189 

.9589 

yv 

2 

.7375 

.8572 

.9254 

.8116 

.8987 

.9474 

N -2 

.6623 

.8083 

.8977 

.7914 

.8861 

.9404 

Table  IX.  Asymptotic  Convergence  Factors  for  the  (2,1)  Cycle  Using  Gauss-Seidel 
Relaxation. 

fact  worth  remembering  when  the  multilevel  solution  is  analyzed. 

Table  IX  indicates  the  asymptotic  convergence  factors  for  a  (2,1)  cycle;  as 
expected,  the  factors  are  lower  than  the  corresonding  factors  for  the  (0, 1)  cycle,  and 
are  universally  lower  than  the  ratios  associated  with  a  single  Gauss-Seidel  sweep  as 
indicated  in  Table  IV.  However,  another  comparison  is  to  raise  the  ratios  in  Table  IV 
to  the  third  power,  corresponding  to  three  relaxation  sweeps  on  the  fine  grid.  In  other 
words,  since  for  a  (2, 1)  cycle  there  are  two  relaxation  sweeps  on  the  fine  grid  before 
the  course  grid  correction,  and  one  afterward,  for  a  total  of  three  relaxation  sweeps 
on  the  fine  grid,  it  seems  reasonable  to  compare  the  performance  of  a  (2, 1)  cycle  to 
three  iterations  of  the  Gauss-Seidel  method,  represented  by  raising  the  Gauss-Seidel 
factor  to  the  third  power.  For  example,  in  Table  IV  for  N  =  32,  j  =  1,  the  ratio 
is  .9567  which,  when  raised  to  the  third  power  is  .8756.  This  value  is  lower  than 
either  of  the  values  in  Table  IX  for  N  =  32,  2j  =  2.  This  also  holds  true  for  the 
asymptotic  convergence  factor  associated  with  N  =  32  and  2j  =  y  (this  phenomenon 
does  not  occur  for  N  =  32,  2j  —  N  —  2,  nor  for  N  =  16).  Thus,  the  overall  predicted 
performance  of  three  Gauss-Seidel  relaxation  sweeps  for  iV  =  32  is  better  than  the 
predicted  performance  of  a  Gauss-Seidel  (2, 1)  cycle. 

Similar  to  the  information  in  Table  VIII,  the  performance  predictions  in  Ta¬ 
ble  IX  improve  with  radial  distance,  and  worsen  with  increase  in  grid  size.  Also,  as 
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(10,10)  cycle,  Gauss-Seidel 

Conservation 

Injection 

N  =  16 

N  =  32 

N  =  64 

II 

N  =  S2 

II 

2 

.3243 

.5716 

.7567 

.3322 

.5751 

.7579 

nv - 

.2678 

.5104 

.7118 

.2948 

.5351 

.7288 

N  —  2 

.2374 

.4795 

.6898 

.2837 

.5256 

.7226 

Table  X.  Asymptotic  Convergence  Factors  for  the  (10,10)  Cycle  Using  Gauss-Seidel 
Relaxation. 

above,  the  conservation  restriction  appears  to  predict  better  performance  than  injec- 
l  ion.  The  factors  are  still  relatively  high,  especially  for  the  N  =  64  grid  size,  so  that 
convergence  will  likely  be  slow. 

The  (10,10)  cycle  is  included  for  reference  so  that  the  performance  of  the 
t.wo-level  cycle  can  analyzed  using  a  large  number  of  iteration  sweeps  on  the  fine  grid. 
The  characteristics  of  the  (10, 10)  cycle  are  similiar  to  those  of  the  (0, 1)  cycle  and  the 
(2, 1 )  cycle;  the  performance  improves  with  radial  distance  and  worsens  with  increase 
in  grid  size.  The  asymptotic  convergence  factors  are  smaller  than  for  previous  cycles, 
I  nit  this  is  expected  since  a  good  deal  more  work  is  represented  by  the  (10, 10)  cycle. 
However,  using  a  comparison  similar  to  that  in  the  discussion  of  Table  IX,  in  Table  IV 
for  N  =  32,  j  =  1,  the  ratio  is  .9567  which,  when  raised  to  the  20th  power  is  .4126. 
Here,  as  in  the  case  of  a  (2,1)  cycle,  this  value  is  lower  than  either  of  the  values 
in  Table  X  for  N  =  32,  j  =  I-  This  holds  true  for  all  the  factors  associated  with 
.V  =  32,  but  is  not  true  for  N  =  16.  Thus,  as  before,  the  predicted  performance 
of  Clauss- .Seidel  relaxation  for  N  =  32  is  better  than  the  predicted  performance  of  a 
Gauss-Seidel  two-level  cycle. 

Tables  XI  and  XII  indicate  the  asymptotic  convergence  factors  for  horizontal 
and  vertical  line  relaxation  (2, 1)  cycles.  Similar  to  the  results  of  the  analysis  on  the 
horizontal  and  line  relaxation  schemes  in  Chapter  IV,  both  line  relaxation  schemes  for 
the  (2, 1)  cycle  indicate  better  performance  than  that  for  the  Gauss-Seidel  method; 
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(2.1)  cycle,  Horizontal  Line  Relaxation 

Conservation 

Injection 

22 

N  =  16 

.V  =  32 

11 

> 

.V  =  16 

II 

N  =  64 

2 

.6784 

.8262 

.9097 

.6938 

.8310 

.9111 

■> 

.6210 

.7836 

.8839 

.6793 

.8203 

.9046 

N  -2 

.5629 

.7408 

.8580 

.6647 

.8096 

.8981 

Ta.l)l('  XI.  .Xsymptotic  Convergence  Factors  for  the  (2,  1)  Cycle  Using  Horizontal  Line 
Rela.vation. 


(2,1)  cycle,  Vertical 

Conservation 

Injection 

2i 

N  =  16 

.V  =  32 

II 

> 

.V  =  16 

iV  =  32 

■V  =  64 

2 

.7258 

.8562 

.9265 

.7429 

.8613 

.9279 

N 

2 

.6326 

.7876 

.8851 

.6928 

.8247 

.9059 

N -2 

.5685 

.7427 

.8586 

.6723 

.8119 

.8988 

Table  XII.  Asymptotic  Convergence  Factors  for  the  (2,1)  Cycle  Using  Vertical  Line 
Relaxation. 

horizontal  line  relaxation  is  the  better  of  the  two.  One  difference  is  that  in  Tables  V 
and  VI,  the  predictions  for  horizontal  line  relaxation  worsen,  and  the  predictions  for 
vertical  line  relaxation  improve,  with  radial  distance,  whereas  in  Tables  XI  and  XII 
the  |)redictions  for  l)oth  line  relaxation  schemes  improve  with  radial  distance.  Using 
calculations  similar  to  those  in  the  discussions  of  Tables  IX  and  X,  the  values  in 
Tables  V  and  VI  are  raised  to  the  third  power  and  compared  with  the  factors  in 
Tables  XI  and  XII.  The  result  is  that  the  predicted  performance  of  three  sweeps  of 
the  horizontal  line  relaxation  is  better  than  that  for  the  (2, 1)  cycle  for  both  N  =  16 
and  /V  =  32.  The  predicted  performance  of  three  sweeps  of  the  vertical  line  relaxation 
is  better  than  that  for  the  (2, 1)  cycle  only  for  N  =  32.  These  results  are  consistent 
with  those  of  the  Gauss-Seidel  relaxation  for  all  of  the  {vi,v'2)  cycles  tested.  That 
is,  for  N  =  32,  the  predicted  performance  of  ordinary  iteration  is  better  than  the 
predictions  for  any  of  the  corresponding  two-level  schemes  considered  here. 
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Ill  the  initial  stages  of  using  V-cycles  to  solve  Equation  IV. 2  at  a  single  time 
step,  we  experiment  with  the  use  of  conservation  restriction  and  linear  interpolation. 
The  result  is  that  this  combination  does  not  result  in  a  coarse  grid  correction.  That 
is.  the  effect  of  applying  the  two-level  scheme  is  to  generate  an  approximation  that 
is  not  as  good  as  that  obtained  by  relaxation  on  the  fine  grid  alone.  Therefore,  we 
experiment  with  the  combination  of  the  injection  restriction  operator  and  the  linear 
interiiolation  operator,  with  the  result  that  this  combination  does  produce  a  coarse 
grid  correction.  .4  correction  scheme  using  V-cycles  with  these  intergrid  transfer 
operators  is  implemented  to  solve  at  a  single  time  step,  and  then  the  solution  is 
stepped  in  time  until  a  “steady”  solution  was  obtained. 

.Multigrid  V-cycle  solutions  to  Equation  IV. 2  are  computed  on  grids  N  =  8, 16, 
and  82.  using  (2.1)  cycles,  Gauss-Seidel  relaxation,  a  =  1  (fully  implicit  time  step- 
])ing),  and  time  step  g  =  /i,  where  the  initial  temperature  distribution  and  boundary 
conditions  are  the  same  as  in  the  iterative  solution.  We  again  use  the  discrete  energy 
norm  (Equation  V.l)  to  measure  the  accuracy  of  these  solutions, 

where  (•,  •)  denotes  the  Euclidean  inner  product,  is  the  operator  defined  in  Equa¬ 
tion  IV. 2,  and  —  T*  approximates  the  discretization  error,  where  is  the 

solution  to  Equation  IV. 2  on  grid  /i,  and  T*  is  the  “exact”  solution  (i.e.,  solution  of 
Equation  IV. 2  with  V  =  64  sampled  on  grid  h)  (see  [Ref.  11]).  The  results  presented 
in  Table  XIII  are  almost  identical  to  those  obtained  for  the  iterative  solution  (see 
Chapter  V). 

As  another  measure  of  how  well  the  iterative  solutions  match  the  multilevel 
solutions,  the  discrete  norm  of  the  solution  differences, 

is  computed  for  each  of  the  above  grid  sizes  and  time  value.  The  results  are  presented 
in  Table  XIV;  while  the  norm  of  the  differences  is  generally  small,  it  increases  with 
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A'  =  8 

A  =  16 

A  =  32 

Initial  time  step 

43.07 

23.97 

6.87 

Time  =  *‘l  second’’  i 
_ 1 

42.99 

23.98 

6.87 

Steady  state 

42.99 

23.98 

6.87 

- 1 

Table  XIII.  Discretization  Errors  for  the  Multilevel  Solution. 


A  =  8 

A  =  16 

A  =  32 

Initial  time  step 

8.00E-16 

1.74E-15 

3.95  E-i  5 

Time  =  "1  second” 

5.68E-16 

1.53E-15 

4.64E-15 

Steady  state 

3.02E-15 

1.02E-14 

5.66E-14 

Table  .XIV.  Norms  of  Differences  B(^tween  Iterative  and  Multilevel  Solutions. 

”  rid  size.  This  might  suggest  a  normalization  problem,  however,  the  factor  of  in  the 
discrete  L'^  norm  is  used  specifically  to  avoid  this  type  of  problem.  More  investigation 
is  required  to  determine  why  there  is  not  better  agreement  between  the  iterative  and 
multilevel  solutions. 

Figures  30  and  31  indicate  the  number  of  iterations  per  time  step  for  the  it¬ 
erative  solutions  and  the  number  of  V-cycles  per  time  step  for  the  multigrid  solution 
for  i\  —  32.  It  is  important  to  note  that,  while  both  processes  require  a  decreasing 
amount  of  work  for  each  time  step  out  to  about  150  steps,  the  multilevel  process  re¬ 
quires  about  I  more  time  steps  than  the  iterative  process  to  reach  a  steady  solution: 
the  multilevel  process  requires  one  V-cycIe  per  time  step  for  time  step  >  150.  In 
addition.  Figures  32  and  33  illustrate  how  the  norm  of  the  difference  between  suc¬ 
cessive  solutions  behaves  as  the  solution  is  stepped  in  time.  As  would  be  hoped,  the 
difference  between  successive  solutions  decreases  with  time  stepping.  However,  as  is 
evident  from  these  latter  figures,  the  time  stepping  process  stalls  somewhat  before 
the  steady  solution  is  reached  (the  tolerance  used  is  machine  e  =  2.22E  —  16).  While 
the  number  of  time  steps  to  reach  a  steady  solution  on  grids  A'  =  8, 16  (not  shown) 
is  about  the  same  for  both  the  iterative  and  multilevel  solutions,  on  grid  N  =  32, 
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as  already  indicated,  the  number  of  time  steps  needed  for  the  multilevel  process  to 
|•('a(•ll  steady  solution  is  about  ^  longer  than  that  required  for  the  iterative  solution. 
Even  so,  the  two  processes  on  all  three  grid  sizes  stall  at  about  the  same  time  step. 


Iterations  per  time  step 


Figure  30.  The  number  of  iterations  per  time  step  for  N  =  32;  170  time  steps  needed 
to  reach  a  steady  solution. 

We  conclude  our  analysis  of  the  multilevel  solution  with  a  discussion  of  how 
nuich  computational  work  must  be  done  to  reach  a  solution.  For  this  purpose,  we 
pick  a  time  step,  for  both  methods,  for  which  the  time  stepping  process  has  stalled. 
For  N  =  8,  use  time  step  s  =  50;  for  N  —  16,  use  time  step  s  =  80;  and  for  N  =  32, 
use  a  time  step  of  s  =  150.  The  norm  of  the  difference  between  successive  solutions 
for  these  time  steps  is  on  the  order  of  10“^®.  We  use  calculations  similar  to  those  in 
[Ref.  3]  in  computing  the  cost.  If  we  consider  a  work  unit,  WU,  to  be  the  cost  of 
one  relaxation  sweep  on  the  finest  grid,  then  we  can  estimate  how  many  work  units 
are  associated  with  each  V-cycle.  Since  we  are  using  a  (2,1)  cycle,  relaxation  occurs 
three  times  on  each  level;  one  relaxation  sweep  on  the  grid  requires  p~‘^  of  a  work 
unit,  where  d  is  the  dimension.  Adding  these  costs,  and  using  the  geometric  series  as 
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V-cycies  per  time  step 


Figure  ;n.  The  number  of  V-cycles  per  time  step  for  N  =  32;  230  time  steps  needed 
to  reach  a  steady  solution. 

an  upper  bound,  vve  have 

V-cycle  Cost  =  3{1  +  2"''  +  •  ■  •  +  2-"''}WU  <  — ^WU. 

Thus,  for  the  two-dimensional  case  we  have  Cost  =  4WU.  A  comparison  of  the 
computational  cost  for  the  various  grid  sizes  is  given  in  Table  XV. 


N  =  S 

N=\Q 

II 

Time  step  (s) 

50 

80 

150 

“Elapsed”  time  (■^) 

5.00 

4.69 

Number  of  V-cycles 

1246 

4076 

13281 

Total  Cost  (WU) 
V-cycle  cost  (WU) 

4984 

16304 

53124 

Iteration  cost  (WU) 

6245 

25631 

103185 

Table  XV.  Computational  Cost  for  the  Iterative  and  Multigrid  Solutions  (Not  Includ¬ 
ing  the  Cost  of  the  Intergrid  Transfers). 
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Norm  of  difference  between  successtve  time  step  solutions 


Figure  32.  Norm  of  difference  between  successive  iterative  solutions  for  N  =  32;  170 
time  steps  needed  to  reach  a  steady  solution. 

The  cost  of  implementing  the  intergrid  transfers  is  not  included  in  the  com- 
l)utations  presented  in  Table  XV.  In  considering  the  amount  of  work  represented 
by  one  work  unit,  counting  “adds”  and  “multiplies”  each  as  one  floating  point  op¬ 
eration  (flop),  the  number  of  flops  for  one  relaxation  sweep  on  a  grid  of  size  N  is 
approximately  41 A^'^.  Injection  restriction  requires  no  arithmetic;  linear  interpola¬ 
tion  requires  approximately  ^  flops.  Thus,  the  cost  of  the  intergrid  transfers  for 
the  two-level  scheme  is  about  {^141N'^)WU  =  ^WU .  Again  using  the  geometric 
series  as  an  upper  bound,  we  have  that  the  cost  of  the  intergrid  transfers  is  given  by 

Transfer  Cost  =  +2-'^  +  ■  ■  ■  +  2-"‘'}WU  <  ^WU. 

328  246 

Thus,  we  update  the  information  in  Table  XV  and  present  the  results  in  Table  XVI. 

While  not  dramatic,  the  use  of  multigrid  does  result  in  a  savings  of  compu¬ 
tational  cost.  The  cost  savings  increase  with  grid  size;  the  savings  on  a  grid  of  size 
=  8  is  only  20%,  on  a  grid  of  size  N  —  32  the  savings  is  almost  50%.  In  consid¬ 
ering  these  results,  several  questions  need  to  be  answered.  For  example,  are  these 
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Norm  of  difference  between  successive  time  step  solutions 


Figure  33.  Norm  of  difference  between  successive  multilevel  solutions  for  N  =  32;  230 
time  steps  needed  to  reach  a  steady  solution. 


N  =  S 

N=16 

N  =  32 

Total  Cost  (WU) 
V-cycle  cost  (WU) 

4992 

16354 

53286 

Iteration  cost  (WU) 

6245 

25631 

Table  XVI.  Computational  Cost  for  the  Iterative  and  Multigrid  Solutions,  Including 
the  C'ost  of  the  Integrid  Transfers. 

results  consistent  with  the  local  mode  analysis  and  the  information  that  comes  from 
the  amplification  matrices?  Does  the  multigrid  solution  show  “typical”  performance; 
does  the  process  work  as  one  would  hope? 

In  general,  the  results  are  consistent  with  the  local  mode  analysis,  which 
predicts  that  convergence  will  be  slow.  The  prediction  of  the  amplification  matri¬ 
ces  is  diflficult  to  apply,  since  the  asymptotic  convergence  factors  apply  to  two-level 
schemes  and  not  to  V-cycles.  Nevertheless,  the  convergence  factors  did  not  predict 
that  the  two-level  scheme  would  provide  improvement  over  the  iterative  process  on 
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grid  size  N  =  32,  and  therefore  these  results  are  somewhat  of  a  pleasant  surprise. 
In  other  words,  multigrid  performed  better  than  might  have  been  expected  based  on 
(lie  two-level  convergence  factors.  However,  typical  multigrid  performance  reaches 
( onvergence  (to  the  level  of  truncation  error,  see  [Ref.  3]  and  [Ref.  2])  in  just  a 
few  V-cycles.  Thus,  the  multilevel  solution  developed  here  is  disappointing,  since 
convergence  requires  a  very  large  number  of  V-cycles. 

There  are  no  obvious  causes  for  the  disappointing  performance  of  the  multi¬ 
level  solution  of  Equation  II. 2.  One  of  the  problems  may  be  the  use  of  cylindrical 
coordinates,  which  results  in  a  radial  bias  in  the  EVE  stencils.  The  discretization  of 
(he  t  ime  derivative  using  finite  differences,  while  using  the  EVE  method  discretize 
the  sjiatial  derivatives,  may  also  have  some  impact  on  the  solution  process.  It  is  not 
clear,  however,  that  either  of  these  affects  the  multilevel  solution  any  differently  than 
it  alfects  the  iterative  solution.  The  choice  of  relaxation  scheme  is  a  major  factor 
in  determining  the  convergence  of  the  process;  perhaps  the  use  of  line  relaxation,  or 
some  other  type  smoother,  will  improve  the  performance.  Again,  however,  it  is  not 
certain  that  this  will  result  in  an  improvement  of  the  multigrid  performance  relative 
to  the  iterative  process.  The  determination  of  the  coarse  grid  and  intergrid  transfer 
operators  is  also  something  of  a  concern.  This  is  one  area  that  clearly  alfects  the 
[performance  of  the  multilevel  solution  as  compared  to  the  iterative  solution.  The  fact 
that  the  predicted  performance  of  iteration  sweeps  alone  is  better  than  the  predicted 
performance  of  the  two-level  schemes  emphasizes  the  need  to  consider  the  use  of  dif¬ 
ferent  coarse  grid  and  intergrid  transfer  operators.  Finally,  the  point-source  problem 
includes  a  singularity  at  the  boundary.  Most  of  the  analysis  has  been  done  for  interior 
points,  and  the  solution  results  do  not  drastically  differ  from  what  is  predicted  by  the 
analysis.  Neverthless,  the  treatment  of  the  boundaries,  particularly  for  this  problem, 
may  provide  a  means  to  improve  the  solution  process. 
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VII. 


CONCLUSION 


The  purpose  of  this  work  is  to  apply  the  FVE  method  to  discretize  the  axisym- 
metric  heat  equation,  and  implement  multigrid  in  solving  it.  In  addition  to  using  the 
FVE  method  to  discretize  the  spatial  derivatives,  finite  differences  are  applied  to  the 
time  derivatives,  resulting  in  a  time-stepping  scheme  that  makes  use  of  a  weighted 
average  of  the  current  and  future  time  steps.  The  use  of  cylindrical  coordinates  results 
in  a  discretization  stencil  that  has  a  radial  bias,  the  effect  of  which  is  evident  in  the 
analysis  of  various  relaxation  schemes.  Gauss-Seidel  is  the  relaxation  scheme  chosen 
ftir  implementation  in  the  iterative  and  multilevel  solution  processes.  Due  to  numer¬ 
ical  anomalies  encountered  in  computing  the  iterative  solution,  the  weighted-average 
time-stepping  scheme  becomes  a  fully  implicit  backward  time-stepping  scheme.  Its 
use  requires  more  computational  work  than  the  weighted-average  scheme,  but  its  use 
is  required  to  achieve  a  physically  meaningful  solution.  The  specifics  of  the  multilevel 
technique  are  then  presented,  including  the  development  of  the  coarse  grid  and  inter¬ 
grid  transfer  operators,  and  the  predictions  of  expected  performance  for  the  two-level 
scheme.  The  multilevel  solution  is  then  implemented;  its  results  are  analyzed  and 
compared  to  the  results  of  the  iterative  solution  process.  The  results  of  our  work  are 
somewhat  disappointing  in  that  we  are  not  yet  able  to  achieve  the  hoped-for  level  of 
success.  There  are,  however,  a  number  of  positive  results  that  come  from  this  work, 
specifically,  the  application  of  the  FVE  method  to  a  problem  in  cylindrical  coordi¬ 
nates  and  the  use  of  Maple  software  to  compute  the  necessary  integrals.  Additionally, 
we  know  that  the  usual  multilevel  techniques  do  not  produce  the  expected  results  and 
therefore  there  is  ample  room  for  further  study. 

There  are  several  possibilities  for  future  work.  One  specific  question  that  has 
not  been  addressed  is  the  determination  of  the  consistency  of  the  discretization  oper¬ 
ator.  Numerical  results  are  useful  in  attempting  to  verify  this  result,  but  consistency 
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slioiild  be  proven  analytically.  In  conjunction  with  this  idea,  there  may  be  a  require¬ 
ment  to  review  the  assumption  that  the  solution  can  be  adequately  represented  by 
a  ri)llection  of  continuous,  piecewise  linear  basis  functions.  There  is  no  evidence  as 
yet  that  would  require  this  review,  but  it  is  possible.  .Additionally,  we  are  assum¬ 
ing  a  uniform  step  size  that  applies  to  both  the  r  and  ^  directions.  One  possible 
modification  to  this  approach  is  to  implement  a  non-uniform  grid,  perhaps  one  that 
results  in  control  volumes  that  are  the  same  for  every  point  on  the  grid.  The  point  of 
t  his  modification  would  be  to  eliminate  the  radial  bias  in  the  discretization.  For  the 
discretization  of  the  time  derivatives,  the  requirement  to  produce  a  physically  mean¬ 
ingful  solution  has  resulted  in  a  fully  implicit  time-stepping  scheme.  While  there  may 
be  another  method  to  discretize  the  time  derivatives,  it  will  have  to  meet  the  same 
criterion  and,  therefore,  may  not  be  much  of  an  improvement. 

All  of  the  above  considerations  apply  generally  to  the  problem,  and  not  specifi¬ 
cally  to  multigrid.  There  are  at  least  three  areas  in  which  multigrid  performance  may 
l>e  improved:  intergrid  transfer  operators,  coarse  grid  operator,  and  treatment  of 
boundary  conditions.  The  intergrid  transfer  operators  used  are  the  simplest  available 
that  result  in  coarse  grid  correction  and,  consequently,  choosing  different  operators 
may  improve  performance.  Other  possibilities  include  the  use  of  a  restriction  operator 
that  is  based  on  the  piecewise  linearity  of  the  basis  functions,  or  perhaps  an  inter- 
|)olation  operator  that  is  based  on  enforcing  conservation.  In  other  words,  choose 
intergrid  operators  that  satisfy  the  variational  condition,  for  c  some 

constant  (see  [Ref.  3]  and  [Ref.  1]).  The  choice  of  the  coarse  grid  operator  may 
need  to  be  made  in  conjunction  with  these  intergrid  operators.  That  is,  despite  the 
computational  complexities  that  would  be  introduced,  choose  a  coarse  grid  operator 
such  that  the  Galerkin  condition  =  (see  [Ref.  3]  and  [Ref.  1])  is  sat¬ 

isfied.  Additionally,  special  treatment  of  the  boundaries  may  be  required,  since  the 
boundary  of  the  point-source  problem  has  a  singularity  at  the  origin  (see  [Ref.  2]). 

Finally,  the  consideration  of  the  point-source  problem  is  the  first  step  in  ap- 
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plying  r,he  FVE  discretization  and  multilevel  solution  to  the  Navier-Stokes  equation 
1  liai  iirises  from  the  molten-pool  problem.  The  intricacies  of  the  molten-pool  problem 
present;  several  challenges.  Included  are  the  requirement  to  solve  a  system  of  three 
PDFs  in  three  unknowns  and  the  requirement  to  keep  track  of  the  moving  phase- 
change  boundary  as  the  molten  pool  expands.  .4dditionally,  high  flow  velocities  and 
small  local  length  scales  result  when  convection  is  vigorous  ([Ref.  13]),  causing  fur- 
1  her  numerical  complications.  It  is  anticipated  that  the  Full  Approximation  .Scheme 
( IfAS)  (see  [Ref.  ‘ij  and  [Ref.  1])  can  be  used  to  effectively  treat  the  non-linear  nature 
of  the  problem,  and  that  the  Fast  .Adaptive  (.Composite  Grid  (FAC)  method  (see  [Ref. 
1]  will  be  useful  in  overcoming  difflculties  that  arise  in  the  geometry  of  the  problem. 
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APPENDIX  A.  THE  DISCRETE  FOURIER 

TRANSFORM 


Below  we  briefly  list  the  definition  and  some  of  the  properties  of  the  Discrete 
Fourier  Transform,  DFT.  For  an  exhaustive  treatment  see  [Ref.  10];  what  appears 
in  this  appendix  closely  follows  that  treatrnent. 

The  DFT  is,  in  at  least  one  interpretation,  a  way  to  discretely  approximate 
the  Fourier  transform.  For  example,  assume  we  are  working  on  a  temporal  or  spatial 
domain  [— 7,  f  j  with  grid  spacing  Ax  =  A/N  and  grid  points  Xj  =  jAx,  and  its 
associated  frequency  domain  [— ^,^]  with  grid  spacing  Au)  =  2t: f A  and  grid  points 
=  in  Au:.  where  in  —  +  1  ;  y .  Suppose  Vj  denotes  the  sampled  values.  n(.r_,  ), 

of  a  function  defined  on  the  domain  [— y,  y],  and  that  is  the  Fourier  transform 

of  V  at  the  frequency  grid  points,  ujm,  where  the  Fourier  transform  is  given  by 

roo 

v{uj)  =  /  v{x)e~'^'^'^^dx. 


Thus,  the  DFT  can  be  considered  to  approximate  the  Fourier  transform  by 

V(u>^)  «  AVm- 

We  now  give  a  formal  definition  of  the  DFTF 


'There  are  several  definitions  of  the  DFT,  as  indicated  in  [Ref.  10]. 
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Definition  A.l  Let  N  be  an  even  positive  integer  and  {vj}  a  sequence'  of  N  complex 
numhers  where  ji  =  —  y  +  1  ;  y.  The  discrete  Fourier  transform  of  the  sequence  {t’j} 
is  another  sequence  of  N  complex  numbers,  {Kn})  whose  elements  are  given  by 


K,.  =  ^  E  » 

-  -T  +  t 


/•  -V  I  I  .  ,v 

Jar  ?;/  =  ——+  1  .  — . 


I'lie  DFT  may  be  defined  for  N  odd  as  well,  however  we  will  not  need  this  definition 
Cor  our  discussion.  The  choice  of  indices,  — -j  +  1  :  y  as  opposed  to  0  :  iV  —  1,  is 
made  deliberately  in  order  to  simplify  some  of  the  analysis  that  is  to  follow  from 
applying  the  DFT,  and  because  it  is  more  natural,  if  less  familiar.  While  the  choice 
of  indices  does  not  affect  the  applicability  of  the  transform,  care  must  be  taken  to 
avoid  confusion  over  which  grid  points  correspond  to  what  type  of  frequency  (see 
Figure  34).  We  use  the  operator  notation  'D{vj}  to  denote  the  DFT  of  the  sequence 
{uj},  and  V{vj}m.  to  indicate  the  mth  element  of  the  transform,  V{vj}m  =  Kn- 

In  general,  the  output  of  the  DFT,  {Kn},  is  a  comp  lex- valued  sequence.  The 
interpretation  of  the  DFT  coefficients  is  essentially  the  same  as  that  given  for  the 
(continuous)  Fourier  transform.  Also,  Kn  is  even  more  closely  related  to  Cm,  the 
Fourier  series  coefficient.  In  many  ways,  the  DFT  is  more  naturally  viewed  as  an 
approximation  to  Cm  than  to  v{uj).  The  mth  DFT  coefficient  Kn  gives  the  “amount” 
of  the  mth  mode  (with  a  frequency  Um)  that  is  present  in  the  input  sequence  Vj.  In 
contrast  to  the  use  of  modes  of  all  frequencies,  as  in  the  Fourier  transform,  an  A-point 
DFT  uses  only  N  distinct  modes,  with  roughly  y  different  frequencies.  The  modes 
can  be  labeled  by  the  frequency  index  m,  and  each  mode  has  a  value  at  each  grid 
point  Xj  where  j  =  — y  +  1  :  y.  Therefore  we  denote  the  jith  component  of  the  mth 


DFT  mode  as 


for  i,m  =  +  1  :  f . 


-jm 

Wjv  =  e  ^  , 


"The  terms  sequence  and  vector  are  used  interchangeably  to  denote  the  input/output  of  the  DFT. 
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Figure  34.  Using  centered  indices  (bottom  figure),  the  low  frequencies  correspond  to 
j/.;|  <  y  for  a  single  set  of  sample  points;  the  high  frequencies  correspond  to  \k\  > 
With  non-centered  indices,  low  frequencies  correspond  to  0  <  ^  and  ^  <  k  <  N 

for  a  single  set;  high  frequencies  correspond  to  ^  <  k  < 

The  DFT  allows  us  to  transform  from  the  temporal  (spatial)  domain  into  the 
frequency  domain.  The  Inverse  Discrete  Fourier  Transform,  IDFT,  permits  us 
to  transform  from  the  frequency  domain  back  into  the  temporal  (spatial)  domain.  We 
now  give  a  definition  for  the  IDFT. 

Definition  A. 2  Let  N  be  an  even  positive  integer  and  {Vm)  ®  sequence  of  N  complex 
numbers  where  m  =  — -y  +  1  :  y.  The  inverse  discrete  Fourier  transform  of  {Vm}  is 
another  sequence  of  N  complex  numbers,  whose  elements  are  given  by 

u,=  t  Kne^,  (A.2) 

m=-f+l 

r  N  i  t  .  N 

.for  J  =  -y  +  1  .  y. 

.‘\s  with  the  DFT,  the  IDFT  is  defined  for  N  odd,  but  we  will  not  need  this  defini¬ 
tion.  Additionally,  we  use  the  operator  notation  to  denote  the  IDFT  of  the 

sequence  {Kn},  and  V~^{Vm}j  to  indicate  the  jth  element  of  the  inverse  transform, 
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The  discrete  Fourier  transform  and  its  inverse  have  many  useful  properties. 
Below  are  listed  some  of  the  more  important  ones  that  we  will  be  using.  First, 
however,  we  give  the  discrete  orthogonality  property  for  the  complex  exponential, 
wliicli  is  essential  in  working  with  DFTs. 


Property  A.l  Orthogonality.  Let  I  and  m  be  integers  and  let  N  be  a  positive 
integer.  Then 


N-l 

Z 


i27nl 

e  ^  e 


t2Trjm 

N 


—  rn), 


(A.3) 


j=0 


inhere  Siwih)  is  the  modular  Kronecker  delta,  defined  bxj 


SN{k)  = 


1  if  k  =  Q  or  a  multiple  of  N 
0  otherwise 


•Although  we  will  not  specifically  need  it,  we  include  the  relationship  between 
the  DFT  and  the  IDFT: 


Property  A. 2  Inversion.  Let  {uj}  be  a  sequence  of  N  complex  numbers  and  let 
{Kn}  the  DFT  of  this  sequence.  Then  {'D{vj}m}j  =  Vj- 

The  remaining  properties  that  we  will  need  are  presented  below.  Proofs  and 
detailed  explanations  are  found  in  [Ref.  10]. 


Property  A.3  Periodicity.  The  complex  sequences  {uj}  and  {Kn}  defined  by  the 
N -point  DFT  pair  (A.l)  and  (A. 2)  are  N -periodic.  That  is, 

Vj  =  Vj±N  and  Vm  =  Kn±iV  for  all  Intergers  j  and  m. 


Property  A. 4  Linearity.  If  {wj}  and  {wj}  are  two  complex-valued  sequences  of 
length  N  and  a  and  j3  are  two  complex  numbers,  then 

V{cxVj  +  =  aV{vj)m  + 


Property  A. 5  Modulation.  If  the  elements  of  the  input  sequence  {vj}  are  multi¬ 
plied  by  for  k  a  fixed  integer,  then  the  DFT  of  the  modulated  sequence  is  given 
by 

=  T>{vj}^.k  =  Vm-k- 
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Property  A. 6  Shifting.  If  the.  input  sequence  {t»j}  is  shifted  (or  translated)  k  U7iits 
to  the  right,  then 

}j7i  — 

VV’e  make  extensive  use  of  the  DFT  in  analyzing  our  solution  process.  The  DFT 
is  an  extremely  useful  tool  in  predicting  the  performance  of  relaxation  schemes  (Chap¬ 
ter  IV)  and  in  evaluating  various  elements  in  the  multilevel  solution  (Chapter  VI)  via 
loccil  mode  analysis.  Of  the  properties  outlined  above,  the  orthogonality  and  shifting 
properties  are  the  two  that  figure  most  prominently  in  local  mode  analysis. 


103 


APPENDIX  B.  INTERPOLATION  TOOLS 


Below  are  the  interpolation  tools  developed  by  David  Canright  [Ref.  8]  for 
(oinputing  the  FVE  stencils  using  Maple.  The  indexing  in  the  following  program 
does  NOT  follow  the  indexing  in  the  text:  the  subscript  i  is  the  column  indicator 
and  corresponds  to  the  column  indicator  j  used  in  the  text;  the  subscript  j  is  the 
row  indicator  and  corresponds  to  the  row  indicator  k  in  the  text.  Thus,  the  indices 
(a.  /)  in  this  Maple  program  indicate  the  jth  row  and  the  ith  column,  in  the  text  the 
indices  {k,j)  indicate  the  A;th  row  and  jth  column. 

The  routine  "plane”  returns  the  function  for  the  plane  through  the  three  given 
point s.  There  is  no  error  checking,  so  it  has  problems  with  special  cases:  all  3  points 
collinear  (infinite  solutions);  and  vertical  planes  (no  solutions). 

plane  :=  proc(xl,yl,zl,x2,y2,z2,x3,y3,z3)  local  a,b,c,x,y; 

unapply( 

subs( 

solve(  { 

a*  x\  +  6*  yl  +  c  =  zl, 
a*  x2-\-b*y2  +  c  =  z2, 
a  *  x'i  +  b  *  y'S  +  c  =  z'i  }, 

{a,  6,  c}), 
a*x  +  b*y  +  c), 

x,y ); 

end: 


To  interpolate  the  indexed  function  z  on  the  6  triangles  surrounding  the  point 
{i  *  h,j  *  fi),  use  the  6  functions  (counter  clockwise  from  northeast)  “tri:  ene,  nne, 
nw,  wsw,  ssw,  se”.  Note:  these  return  functions  of  two  variables  (for  example  r,z). 
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triene  :=  (i,j,z)  -  >  planefi=*'h.j^li.z._p,(i+l)*h.j*h,z.-e,(i+l)*h, 

(j  +  i  )*h,z.-iie): 

trinne  :=  (i,j,z)  -  >  plane(i=*'h,j*h,z._p,i^h,(j+l)*h,z._n,(i+l)*li, 
(j  +  l)=^h,z.-iie): 

trinw  :=  (i,j,z)  -  >  plane(i^hj="h.z._p,i="h.(j+l)*h,z._n,(i-l)*h, 

triwsw  ;=  (i,j,z)  -  >  plane(i*h,j’"h,z._p,(i-l)*li,j*h,z.-w,(i-l)*h, 
(j-l)*h,z._sw): 

trissw  ;=  (i,j,z)  —  >  plane(i*h,j'*'h,z._p,i*h,(j  — l)*h,z.-s,(i— l)*h, 
(j-l)*h,z._sw): 

trise  :=  (i,j,z)  -  >  plane(i*h,j'*'h,z._p,i*h,(j-l)’^h,z.j,(i+l)'*'h, 
j^h,z._e): 


To  convert  the  notation  from  subscripted  back  to  indexed,  use  "toindex(expression. 
variable)”,  where  the  variable  is  to  become  indexed. 


toindex  :=  proc  (expr,var)  local  n,m; 
subs( 

zip(  (x,y)  -  >  (x=y), 

[var.  (_sw,_s,_se,_w,_p,_e,_nw,_Q,_ne)], 

[seq  (seq  (seq  (var[i+n,j+m],  n=  — l..l),m=  — 1..1)]), 

expr  ); 
end: 


Now  we  define  shorthand  for  limits  of  integration  about  the  general  point 
Here,  rdiag  means  the  diagonal  line,  where  r  depends  on  2.  (This  assumes  integrating 
hrst  in  r,  then  in  z.) 

rp  :=  i*h;  rw  :=  (i-l/2)*h;  re  :=  (i+l/2)*h;  rdiag  :=  i*h-(z-j*h); 
zp  :=  j*h;  zs  :=  (j-l/2)’^h;  zn  :=  (j+l/2)*h; 


Now  compute  the  volume  integral  for  the  unknown  temperature  T  (which  is  a 
surface  integral  of  T  after  factoring  out  the  27r  from  the  ip  integration): 

int  (int  (triene  (i  j,T)(r,z)*r,  r=rdiag..re),z=zp..zn)  + 
int  (int  (trinne  (io,T)(r,z)*r,  r=rp..rdiag),z=zp..zn)  + 
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int  (int  (trinw  (i,j,T)(r,z)*r,  r=rw..rp),z=zp..zn)  + 
int  (int  (triwsvv  (i,j,T)(r,z)*r,  r=rw..rdiag),z=zs..zp)  + 
int  (int  (trissw  (i,j,T)(r,z)'"r,  i'=rcliag..rp).z=zsi.zp)  + 
int  (int  (trise  (i,j,T)(r,z)''r.  r=rp..re),z=zs..zp); 

factor(simplify(” )); 
toindex(”,T); 


'Phis  set  of  computations  produces  the  stencil  entries  for  the  volume  integral: 

^((;P2i>  +(32z>5)r,^,_,  +224rn,  +(32f-  ii)r^ 

(16i  —  +  (32i  —  5)Tjj+i  +  (16i  4-  5)T,+i,j+i) . 

In  similar  fashion,  the  surface  integral  of  VT  is  computed  (which  is  a  circula¬ 
tion  integral  of  T  after  factoring  out  the  27r  from  the  (p  integration): 

int  (+D[1]  (triene(i,j,T))(re,z)*re,z=zp..zn)  -f 
int  (4-D[2]  (trinne(i,j,T))(r,zn)*r,r=rp..re)  + 
int  (-1-D[2]  (trinw(i,j,T))(r,zn)*r,r=rw..rp)  -f 
int  (-D[l]  (trinw(i,j,T))(rw,z)*rw,r=zp..zn)  + 
int  (-D[l]  (triwsw(i,j,T))(rw,z)*rw,z=zs..zp)  -t- 
int  (-D[2]  (trissw(i,j,T))(r,zs)*r,r=rw..rp)  -f 
int  ^D[2]  (trise(i,j,T))(r,zs)=^r,r=rp..re)  4- 
int  (4-D[l]  (trise(i,j,T))(re,z)*re,z=zs..zp); 

factor(simplify(” )); 
toindex(”,T); 


This  set  of  computations  produces  the  stencil  entries  for  the  surface  integral: 


h 

-  ((1  4-  20r.+, j  -  8tT,j  -h  (2i  -  +  2tT,j+i  + 
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